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Abstract. In this paper we present efficient algorithms for the com- 
putation of several invariant objects for Hamiltonian dynamics. More 
precisely, we consider KAM tori (i.e diffeomorphic copies of the torus 
such that the motion on them is conjugated to a rigid rotation) both La- 
grangian tori (of maximal dimension) and whiskered tori (i.e. tori with 
hyperbolic directions which, together with the tangents to the torus and 
the symplectic conjugates span the whole tangent space). In the case of 
whiskered tori, we also present algorithms to compute the invariant split- 
ting and the invariant manifolds associated to the splitting. We present 
them both for the case of discrete time and for differential equations. 

The algorithms for tori are based on a Newton method to solve 
an appropriately chosen functional equation that expresses invariancc. 
Among their features we highlight: 

• The algorithms are efficient: if we discretize the objects by N ele- 
ments, one step of the Newton method requires only O(N) storage 
and 0(N ln(iV)) operations. Furthermore, if the object we con- 
sider is of dimension £, we only need to compute functions of I 
variables, independently of what is the dimension of the phase 
space. 

• The algorithms do not require that the system is presented in 
action-angle variables nor that it is close to integrable. 

• The algorithms are backed up by rigorous a-posteriori bounds 
which state that if the equations are solved with a small resid- 
ual and some explicitly computable condition numbers are not too 
big, then, there is a true solution which is close to the computed 
one. 

• The algorithms apply both to primary (i.e non-contractible) and 
secondary tori (i.e. contractible to a torus of lower dimension, such 
as islands). They also apply to whiskered tori. 

The algorithms for invariant splittings are based on computations of 
proyections (rather than in graph transforms). The computations of 
invariant manifolds are also efficient in the sense indicated before. 

The algorithms we present have already been implemented. We will 
report on the technicalities of the implementation and the results of 
running them elsewhere. 
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1. Introduction 

The goal of this paper is to present efficient algorithms to compute very ac- 
curately several objects of interest in Hamiltonian dynamical systems (both 
discrete-time dynamical systems and differential equations) . More precisely, 
we present algorithms to compute: 

• Lagrangian KAM tori. 

• Whiskered KAM tori. 

• The invariant bundles of the whiskered tori. 

• The stable and unstable manifolds of the whiskered tori. 

The algorithms are very different. For example, the algorithms for tori 
require the use of small divisors and symplectic geometry and the algorithms 
for invariant bundles and invariant manifolds rely on the theory of normal 
hyperbolicity and dichotomies. The computation of whikered tori has to 
combine both. 

We recall that KAM tori are manifolds diffeomorphic to a torus which 
are invariant for a map or flow, on which the motion of the system is con- 
jugate to a rotation. As we will see later, this is also equivalent to quasi- 
periodic solutions. The tori are called Lagrangian when they are Lagrangian 
manifolds, which in our case amounts just to the fact that the tori have a 
dimension equal to the number of degrees of freedom of the system. The 
tori are called whiskered when the linearized equation has directions that 
decrease exponentially either in the future (stable) or in the past (unstable) 
and these directions together with the tangent to the torus and its symplec- 
tic conjugate span the whole tangent space. These invariant spaces for the 
linearization have non-linear analogues, namely invariant manifolds. It has 
been recognized since |Arn64| that whiskered tori and their invariant man- 
ifolds are very interesting landmarks that organize the long-term behavior 
of many systems. 

The algorithms we present are based on running an efficient Newton 
method to solve a functional equation, which expresses the dynamical prop- 
erties above. What we mean by efficient is that if we discretize the problem 
using N Fourier coefficients, we require 0(N) storage and only 0(Nln(N)) 
operations for the Newton step. Since the functions we are considering are 
analytic, we see that the truncation error is 0(exp(— CN 1 ^)) where d is 
the dimension of the object. Note that, in contrast, a straightforward im- 
plementation of a Newton method would require to use 0(N 2 ) storage - 
to store the linearization matrix and its inverse - and 0(N 3 ) operations to 
invert. 

In practical applications, using the algorithms described in this paper, 
computing with several million coefficients becomes quite practical in a typ- 
ical desktop computer of today. Implementation details and the results of 
several runs will be discussed in another companion paper [HdlLS09] . Given 
the characteristics of today's computers, savings in storage space are more 
crucial than savings in operations for these problems. 
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The algorithms we present here are inspired by the rigorous results of 
[LG.TV05] - for KAM tori - and [FdlLS09bl IFdlLSOflaj -for whiskered tori. 
The algorithms to compute the stable and unstable manifolds had not been 
previously discussed. The rigorous results of the above papers are also based 
on a Newton method applied to the same functional equation that we con- 
sider here. Of course, going from a mathematical treatment to a practical 
algorithm requires significantly many more details. In particular, the algo- 
rithms to compute invariant splittings and invariant manifolds are different 
from those in the above references. This paper discusses these algorithmic 
issues. 

The results of the papers [LGJV051 IFdlLS09a] give a justification of 
the algorithms for tori and splittings presented here. The theorems in 
[LGJV05, FdlLS09a], have been formulated in an a-posteriori way, i.e. the 
theorems assert that if we have a function which solves approximately the 
invariance equation very accurately (e.g. the outcome of a successful run of 
the algorithms) and which also satisfies some explicit non-degeneracy con- 
ditions, then, we can conclude that there is a true solution which is close to 
the computed solution. Hence, by supplementing our calculations with the 
(very simple) computations of the non-degeneracy conditions (they play a 
role very similar to the condition numbers common in numerical analysis), 
we can be sure that the computation that we are performing is meaningful. 
This allows to compute with confidence even close to the limit of validity of 
the KAM theorem (a rather delicate boundary since the smooth KAM tori 
do not disappear completely but rather morph into Cantor sets). 

Since the papers [LGJV05} FdlLS09a] contain estimates, in the present 
paper, we will only discuss the algorithmic issues. For example, we will detail 
how solutions of equations (whose existence was shown in the above papers) 
can be computed with small requirements of storage and small operation 
count. Note that different algorithms of a same mathematical operation can 
have widely different operation counts and storage requirements. (See, for 
example, the discussion in |Knu97| on the different algorithms to multiply 
matrices, polynomials, etc.) On the other hand, we will not include some 
implementation issues (methods of storage of arrays, ordering of loops, pre- 
cision, etc.) needed to obtain actual results in a real computer. They will 
be given in another paper together with experimental results obtained by 
running the algorithms. 

One remarkable feature of the algorithms presented here is that they do 
not require the system to be close to integrable. We only need a good 
initial guess for the Newton method. Typically, one uses a continuation 
method starting from an integrable case, where solutions can be computed 
analytically. However, in the case of secondary KAM tori, which do not exist 
in the integrable case, one can use, for instance, Lindstedt series, variational 
methods or approximation by periodic orbits to obtain an initial guess. 
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As for the algorithms to compute invariant splittings, we depart from the 
standard mathematical methods (most of the time based on graph trans- 
forms) and we have found more efficient to device an equation for the invari- 
ant projections. We also found some acceleration of convergence methods 
that give superexponential convergence. They are based on fast algorithms 
to solve cohomology equations which could be of independent interest (See 
Appendix lAl). 

The algorithms to compute invariant manifolds are based on the parame- 
terization method |CFL03a| ICFL05| . Compared to standard methods such 
as the graph transform it has the advantage that to compute geometric ob- 
jects of dimension £, we only need to compute with functions of dimension 
1. In contrast with |CFL03a] ICFL05] . which was based on contractive itera- 
tions, our method is based on a Newton iteration which we also implement 
without requiring large matrices and requiring only Nlog(N) operations. 

An overview of the method. The numerical method we use is based 
on the parameterization methods introduced in [CFL03a, CFL03b]. In this 
section, we provide a sketch of the issues, postponing some important details. 
For brevity, we make the presentation for maps only, even if a similar sketch 
can be made for flows. 

Invariant tori. We observe that if F is a map and we can find an embedding 
K in which the motion on the torus is a rotation u, it should satisfy the 
equation 

F(K(9))-K(0 + u) = 0. (1) 
Given an approximate solution of J!]), i.e. 

F(K(9))-K(0 + uj) = E(0), 

the Newton method aims to find A solving 

DF{K(6))A{6)-A(6 + lo) = -E{8), (2) 

so that K + A will be a much more approximate solution. 

The main idea of the Newton method is that, using the decomposition 
into invariant subspaces, one can decompose ([2]) into three components 

DF(K(8))A S (8) - A s (6 + u) = -E s (6) 

DF(K{e))A u {e)-A u {8 + oo) = -E u {9) (3) 

DF(K(0))A C (8) - A c (8 + u) = -E c (6) 

where the s, u refer to the stable, unstable components and c refers to the 
component along the tangent to the torus and its symplectic conjugate. For 
Lagrangian tori, only the E c part appears in the equations. 
The algorithm requires: 

• Efficient methods to evaluate the LHS of ([I]). 

• Efficient methods to compute the splitting. 

• Efficient methods to solve the equations ([3]). 
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As we will see in Section [3T6l to evaluate ((TJ), it is efficient to use both a 
Fourier representation (which makes easy to evaluate K{6 + w)) and a real 
representation which makes easy to evaluate F(K(9)). Of course, both of 
them are linked through the Fast Fourier Transform (FFT from now on). 

The methods to compute the splitting are discussed in Section 14.31 More 
precisely, we present a numerical procedure to compute the projections 
on the linear stable/unstables subspaces based on a Newton method. In 
[HdlLS09], we present an alternative procedure for the computation of the 
projections based on the calculation of invariant bundles for cocycles. In- 
deed, these algorithms require the computation of the projections over the 
linear subspaces of the linear cocycle. 

The solution of the hyperbolic components in equation ([3]) is discussed 
in Section 14.11 and Appendix [A] Indeed, equations of this form appear as 
well in the calculation of the invariant splitting discussed in Section 14.31 A 
first method is based on an acceleration of the fixed point iteration (Ap- 
pendix [AjJ. We note that to obtain super exponential convergence for the 
solution of ([3]), we need to use both the Fourier representation and the real 
space representation. 

In the case that the bundles are one-dimensional, there is yet another 
algorithm, which is even faster than the previous ones (see Appendix IA.2|) , 
The algorithms are discussed for maps, and they do not have an easy analog 
for flows except by passing for the integration of differential equations. We 
think that this is one case where working with time-1 maps is advantageous. 

The most challenging step is the solution of the center component of ([3j). 
This depends on cancelations which use the symplectic structure, involves 
small divisors and requires that certain obstructions vanish. Using several 
geometric identities that take advantage of the fact that the map is symplec- 
tic (see Section I4.2p , the solution of ([3]) in the center direction is reduced to 
solving the following equation for ip given rj, 

p(6)-p(e + u)=rj(e). (4) 

Equation (j3J) can be readily solved using Fourier coefficients provided that 
J r] = (and that uj is sufficiently irrational). The solution is unique up to 
addition of a constant. 

The existence of obstructions - which are finite dimensional - is one of 
the main complications of the problem. It is possible to show that, when the 
map F is exact symplectic, the obstructions for the solution are 0(||£?|| 2 ). 
An alternative method to deal with these obstructions is to add some new 
- finite dimensional - unknowns A and solve, instead of ([I]), the equation 

F(K(6) - K(9 + + G(9)X = 

where G{9) is an explicit function. Even if A is kept through the iteration 
involving approximate solutions, it can be shown that, if the map is exact 
symplectic, we have A = 0. This counterterm approach also helps to weaken 
non-degeneracy assumptions. 
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A minor issue that we omit in this preliminary discussion is that the 
solutions of ([I]) are not unique. If K is a solution, K defined by K(9) = 
K{6 + a) is also a solution for any a G Mr. This can be easily solved by 
taking an appropriate normalization that fixes the origin of coordinates in 
the torus. In [FdlLS09a] it is shown that this is the only non-uniqueness 
phenomenon of the equation. Furthermore, this local uniqueness property 
allows to deduce results for vector fields from the results for maps. 

It is important to remark that the algorithms that we will present can 
compute in a unified way both primary and secondary tori. We recall here 
that secondary tori are invariant tori which are contractible to a torus of 
lower dimension, whereas this is not the case for primary tori. The tori which 
appear in integrable systems in action-angle variables are always primary. In 
quasi- integrable systems, the tori which appear through Lindstedt series or 
other perturbative expansions starting from those of the integrable system 
are always primary. Secondary tori, however, are generated by resonances. 
In numerical explorations, secondary tori are very prominent features that 
have been called "islands". In [HLOOJ, one can find arguments showing that 
these solutions are very abundant in systems of coupled oscillators. As an 
example of the importance of secondary tori we will mention that in the 
recent paper [DLS06] they constituted the essential object to overcome the 
"large gap problem" and prove the existence of diffusion. In [DH09], one 
can find a detailed analysis of secondary tori. 

In this paper, we will mainly discuss algorithms for systems with dynamics 
described by diffeomorphisms. For systems described through vector fields, 
we note that, taking time— 1 maps or, more efficiently, surfaces of section, 
we can reduce the problem with vector fields to a problem with diffeomor- 
phisms. However, in some practical applications, it may be convenient to 
have a direct treatment of the system described by vector fields. For this 
reason, we have included the invariance equations for flows, in parallel with 
the invariance equations for maps and we have left for the Appendix the 
algorithms that are specially designed for flows. 

Invariant manifolds attached to invariant tori. When the torus is whiskered, 
it has invariant manifolds attached to it. For simplicity, in this presentation 
we will discuss the case of one dimensional directions - even if the torus can 
be of higher dimension. 

We use again a parameterization method. Consider an embedding W in 
which the motion on the torus is a rotation oj and the motion on the stable 
(unstable) whisker consists of a contraction (expansion) at rate fx, it should 
satisfy the invariance equation 

F(W(6,s))-W{9 + uj,fis) = 0. (5) 

Again, the key point is that taking advantage of the geometry of the 
problem we can devise algorithms which implement a Newton step to solve 
equation ([5]) without having to store — and much less invert — a large matrix. 
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We first discuss the so-called order by order method, which serves as a 
comparison with more efficient methods based on the reducibility. Although 
they are based on the same idea as for the case of tori, they have not been 
introduced previously and constitute one of the main novelties of this paper. 
We present algorithms that compute at the same time the torus and the 
whiskers and algorithms that given a torus and the linear space compute 
the invariant manifold tangent to it. It is clearly possible to extend the 
method to compute stable and unstable manifolds in general dimensions (or 
even non-resonant bundles). To avoid increasing the length of this paper 
and since higher dimensional examples are harder numerically, we postpone 
this to a future paper. 

Some remarks on the literature. Invariant tori in Hamiltonian dynam- 
ics have been recognized as important landmarks in Hamiltonian dynamics. 
In the case of whiskered tori, their manifolds have also been crucial for the 
study of Arnold diffusion. 

Since the mathematical literature is so vast, we cannot hope to summarize 
it here. We refer to the rather extensive references of [LlaOl] for Lagrangian 
tori and those of [FdlLS09a] for whiskered tori. We will just briefly mention 
that [Gra741 IZeh76] the earliest references on whiskered theory, as well as 
most of the later references, are based on transformation theory, that is 
making changes of variables that reduce the perturbed Hamiltonian to a 
simple form which obviously presents the invariant torus. From the point of 
view of numerics, this has the disadvantage that transformations are very 
hard to implement. 

The numerical literature is not as broad as the rigorous one, but it is still 
quite extensive. The papers closest to our problems are |HL06cI [HL06b, 
HL07J, which consider tori of systems under quasi-periodic perturbation. 
These papers also contain a rather wide bibliography on papers devoted to 
numerical computation of invariant circles. Among the papers not included 
in the references of the papers above because these appeared later, we men- 
tion |CdlL09j . which presents other algorithms which apply to variational 
problems (even if they do not have a Hamiltonian interpretation) . Another 
fast method is the "fractional iteration method" [SimOO]. Note that the 
problems considerd in jHL06cl iHL06bl IHL07] do not involve center direc- 
tions (and hence, do not deal with small divisors) and that the frequency 
and one of the coordinates of the torus is given by the external perturba- 
tion. The methods of [HL06c, HL06b, HL07] work even if the system is not 
symplectic (even if they can take advantage of the symplectic structure). 

The papers [JO051 IJO09] present and implement calculations of reducible 
tori. This includes tori with normally elliptic directions. The use of re- 
ducibility indeed leads to very fast Newton steps, but it still requires the 
storage of a large matrix for the changes of variables. As seen in the exam- 
ples in [HL071 lHL06aj . reducibility may fail in a codimension 1 set (a Cantor 
set of codimension 1 manifolds for elliptic tori in Hamiltonian systems). For 
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these reasons, we will not discuss methods based on reducibility in this pa- 
per (even if it is a useful and practical tool) and just refer to the references 
just mentioned. Indeed, thanks to hyperbolicity, reducibility is not needed 
in the present paper. 

The paper is organized as follows. In Section [2] we summarize the notions 
of mechanics and symplectic geometry we will use. In Section[3]we formulate 
the invariance equations for the objects of interest (invariant tori, invariant 
bundles and invariant manifolds) and we will present some generalities about 
the numerical algorithms. 

Algorithms for whiskered tori are discussed in Section HI In particular, 
we discuss how to compute the decomposition ([3|) of the linearized equation 
(|2|) , and how to solve efficiently each equation in ([3]) . 

In Section [5] we discuss fast algorithms to compute rank-1 (un)stable 
manifolds of whiskered tori. More precisely, we present an efficient Newton 
method to solve equation ©. 

In Appendix [A] one can find the fast algorithms to solve cohomology 
equations with non-constant coefficients that will be used in the computation 
of the splitting ([3]) as well as to solve the hyperbolic components of equations 
0. In Appendices iBllEl one can find the algorithms specially designed for 
flows, analogous to the ones for maps. 

2. Setup and conventions 

We will be working with systems defined on an Euclidean phase space 
endowed with a symplectic structure. The phase space under consideration 
will be 

M C R 2d ' e x T e . 

We do not assume that the coordinates in the phase space are action-angle 
variables. Indeed, there are several systems (even quasi-integrable ones) 
which are very smooth in Cartesian coordinates but less smooth in action- 
angle variables (e.g., neighborhoods of elliptic fixed points [FGB98, GFB98J, 
hydrogen atoms in crossed electric and magnetic fields [RC95, R C97] and 
several problems in celestial mechanics [CC07| ). 

We will assume that the Euclidean manifold Ai is endowed with an exact 
symplectic structure VL = da (for some one-form a) and we have 

U z (u,v) = {u,J(z)v), 

where (•, ■} denotes the inner product on the tangent space of M. and J(z) 
is a skew-symmetric matrix. 

An important particular case is when J induces an almost-complex struc- 
ture, i.e. 

J 2 = — Id . (6) 

Most of our calculations do not need this assumption. One important case, 
where the identity ([6]) is not satisfied, is when J is a symplectic structure on 
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surfaces of section chosen arbitrarily in the energy surface or when J is the 
symplectic form expressed in symplectic polar coordinates near an elliptic 
fixed point. When ([6]) holds, some calculations can be made faster. 

As previously mentioned, we will be considering systems described either 
by diffeomorphisms or by vector-fields. 

2.1. Systems described by diffeomorphisms. We will consider maps 
F :U C M i—t- Ai which are not only symplectic (i.e. F*0, = f2) but exact 
symplectic, that is 

F*a = a + dP, 

for some smooth function P, called the primitive function. 

We will also need Diophantine properties for the frequencies of the torus. 
For the case of maps, the useful notion of a Diophantine frequency is: 

V(u, t) = iu G R £ | \u) ■ k - n^ 1 < v\k\ T V k G Z e - {0}, n G z| , v > L 

2.2. Systems described by vector fields. We will assume that the sys- 
tem is described by a globally Hamiltonian vector-field X, that is 

X = JVH, 

where H is a globally defined function on T*A4. 

In the case of flows, the appropriate notion of Diophantine numbers is: 

£> aff (V, r ) = j w €H* | |w ■ k\- x < v\k\ r V k G J} - {0}} , v>l-\ 

Remark 1. It is well known that for non-Diophantine frequencies substan- 
tially complicated behavior can appear jl ler!)2, iFKWOl] ■ Observing convinc- 
ingly Liouvillian behaviors seems a very ambitious challenge for numerical 
exploration. 

3. Equations for invariance 

In this section, we discuss the functional equations for the objects of inter- 
est, that is, the invariant tori and the associated whiskers. These functional 
equations, which describe the invariance of the objects under consideration, 
are the cornerstone of the algorithms. 

3.1. Functional equations for whiskered invariant tori for diffeo- 
morphisms. At least at the formal level, it is natural to search quasi- 
periodic solutions with frequency uj (independent over the integers) under 
the form of Fourier series 

X M = ^ x k e^ ik ^ n , (7) 
kez 1 

where w G M f and n G Z. 

We allow some components of x in (|7|) to be angles. In that case, it suffices 
to take some of the components of x modulo 1. 
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It is then natural to describe a quasi-periodic function using the so-called 
"hull" function K : T e -> M defined by 

K{6) = J2 x k e 2mk - e , 

so that we can write 

x< n > = K{nuj). 

The geometric interpretation of the hull function is that it gives an em- 
bedding from into the phase space. In our applications, the embedding 
will actually be an immersion. 

It is clear that quasi-periodic functions will be orbits for a map F if and 
only if the hull function K satisfies: 

FoK-KoT lv = 0, (8) 

where denotes a rigid rotation 

T u (9) = 9 + u. (9) 

A modification of the invariance equations ([8]) which we will be important 
for our purpose consists in considering 

FoK-KoT^- (JiKo^DKo) o T W X = 0, (10) 

where the unknowns are now K : T e — > Ai (as before) and A € Mr. Here, Kq 
denotes a given approximate (in a suitable sense which will be given below) 
solution of the equation (JSj) . 

It has been shown in [FdlLS09b, FdlLS09aJ (the vanishing lemma that, 
for exact symplectic maps, if {K, A) satisfy the equation (|10p with Kq close 
to K, then at the end of the iteration of the Newton method, we have 
A = and, therefore, K is a solution of the invariance equation ([8]). In 
other words, the formulations (|10p and ([8]) are equivalent. Of course, for 
approximate solutions of the invariance equation ([8]) , there is no reason why 
A should vanish and it is numerically advantageous to solve the equation 
with more variables. 

The advantage of equation (I10p for numerical calculations is that, at the 
initial stages of the method, when the error in the invariance equation is 
large, it is not easy to ensure that certain compatibility conditions (see (j24"j) 
in Section I3.T[) . are satisfied approximately, so that the standard Newton 
method has problems proceeding. On the other hand, we can always proceed 
by adjusting the A. This is particularly important for the case of secondary 
tori that we will discuss in Section 13. 41 We also note that this procedure 
makes possible to deal with tori when the twist condition degenerates. 

The equations ([8j) and (jlOh will be the centerpiece of our treatment. We 
will discretize them using Fourier series and study numerical methods to 
solve the discretized equations. 

It is important to remark that there are a posteriori theorems (see [FdlLS09b, 
FdlLS09a]) for equations ([8]), (fTUl) (as well as their analogous for flows (fTTI) . 
(|13p ). That is, theorems that ensure that given a function that satisfies 
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(JH|), (|10p up to a small error and that, at the same time, satisfies some non- 
degeneracy conditions (which are given quite explicitly) , then there is a true 
solution close to the computed one. Hence, if we monitor the non- degeneracy 
conditions, we can be sure that the computed solutions correspond to some 
real effects and are not spurious solutions. 

Remark 2. Notice that for whiskered tori the dimension of the torus I is 
smaller than half the dimension of the phase space 2d. Hence, the algorithms 
presented here have the advantage that they look for a function K of £ 
variables to compute invariant objects of dimension t. This is important 
because the cost of handling functions grows exponentially fast with the 
number of variables. Indeed, to discretize a function of i variables in a grid 
of side h into R M , one needs to store (l//i)^ • 2d real values. 

Remark 3. Equations ([8]) and (fTUj) do not have unique solutions. Observe 
that if K is a solution, for any a € R , K o T a is also a solution. In 
[FdlLS09a] , it is shown that, in many circumstances, this is the only non 
uniqueness phenomenon in a sufficiently small neighborhood of K. Hence, 
it is easy to get rid of it by imposing some normalization. See Section 13.5.21 
for a discussion on this issue. 

3.2. Functional equations for whiskered invariant tori for vector- 
fields. In this case, one can write 

k& 1 

where u € R , t € R and then the hull function K is defined by 

x(t) = K(ut). 
The invariance equation for flows is written: 

d^K -XoK = 0, (11) 

where d u denotes the derivative in direction ui 

e 

doj = ^2^kde k - (12) 

k=l 

The modification of (llip incorporating a counterterm is: 

d^K-XoK- JiKo^iDX o K )X = 0, (13) 

where Kq is a given embedding satisfying some non-degeneracy conditions. 

Remark 4. Recall that, taking time— 1 maps, one can reduce the problem 
of vector fields to the problem of diffeomorphisms. Furthermore, since au- 
tonomous Hamiltonian systems preserve energy, we can take a surface of 
section and deal with the return map. This reduces by 2 the dimension of 
the phase space and the parameterization of the torus requires 1 variable 
less. In practice, it is much more efficient to use a numerical integrator to 
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compute the point of intersection with the surface of section than to deal 
with functions of one more variable and with two more components. 

3.3. Some global topological considerations. In our context, both the 
domain and the range of K have topology. As a consequence, there 
will be some topological considerations in the way that the torus T £ gets 
embedded in the phase space. More explicitly, the angle variables of T e can 
get wrapped around in different ways in the phase space. 

A concise way of characterizing the topology of the embedding is to con- 
sider the lift of K to the universal cover, i.e. 

K : R e -+ R 2d ~ e x R e , 

in such a way that K is obtained from K by identifying variables in the 
domain and in the range that differ by an integer. 
It is therefore clear that VeeZ 1 

K p (9 + e) = K p (9), 

( 14 ) 

K q (6 + e) = K q (e) + I(e), 

where K p ,K q denote the projections of the lift on the p and q coordinates 
of M. 2d - £ xl'. It is easy to see that /(e) is a linear function of e, namely 

J(e)i=W= feW) (15) 
\ j=1 / i=i,...,e 

with Iij € Z. 

We note that if a function K q satisfies 

K q (9 + e) = K q (9) + 1(e) , 

the function 

K q (9) = K q (9) - 1(9) (16) 
is e— periodic. The numerical methods will always be based on studying the 
periodic functions K q , but we will not emphasize this unless it can lead to 
confusion. 

Of course, the integer valued matrix I = {Iij}^ remains constant if we 
modify the embedding slightly. Hence, it remains constant under continuous 
deformation. For example, in the integrable case with £ = d, invariant tori 
satisfy K q (9) = 9, so that we have I = Id. Hence, all the invariant tori which 
can be continued from tori of the integrable system will also have I = Id. 

3.4. Secondary tori. One can produce other ^-dimensional tori for which 
the range of / is of dimension less than I. These tori are known as secondary 
tori. It is easy to see that if rank(J) < £ we can contract K(T l ) to a 
diffeomorphic copy of T rank (^ . Even in the case of maximal tori £ = d, 
one can have contractible directions. The most famous example of this 
phenomenon are the "islands" generated in twist maps around resonances. 
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Secondary tori do not exist in the integrable system and they cannot be 
even continuously deformed into some of the tori presented in the integrable 
system. This is often described informally as saying that the secondary tori 
are generated by the resonances. 

Perturbative proofs of existence of secondary tori are done in [LW04J and 
in [DLS06] and in more detail in [DH09]. In [Dua94j one can find rigorous 
results showing that these islands have to be rather abundant (in different 
precise meanings) in many classes of 2D-maps. In particular, for standard- 
like maps, secondary tori appear at arbitrarily large values of the parameter. 

In [HLOO], there are heuristic arguments and numerical simulations argu- 
ing that in systems of coupled oscillators, the tori with contractible directions 
are much more abundant than the invariant tori which can be continued from 
the integrable limit. 

In view of these reasons, we will pay special attention to the computation 
of secondary tori. 

We want to emphasize on some features of the method presented here, 
which are crucial for the computation of secondary tori: 

• The method does not require neither the system to be close to inte- 
grable nor to be written in action-angle variables. 

• The modification of the invariance equations (|8|) and allows to 
adjust some global averages required to solve the Newton equations 
(see |FdlLS09a] ). 

• The periodicity of the function K can be adjusted by the matrix I 
introduced in (|14p . Hence, the rank of the matrix I has to be chosen 
according to the number of contractible directions. 

3.5. Equations for the invariant whiskers. Invariant tori with I < d 
may have associated invariant bundles and whiskers. We are interested in 
computing the invariant manifolds which contain the torus and are tangent 
to the invariant bundles of the linearization around the torus. This includes 
the stable and unstable manifolds but also invariant manifolds associated 
to other invariant bundles of the linearization, such as the slow manifolds, 
associated to the less contracting directions. 

Using the parameterization method, it is natural to develop algorithms 
for invariant manifolds tangent to invariant sub-bundles that satisfy a non- 
resonance condition (see [CFL03a]). This includes as particular cases, the 
stable/unstable manifolds, the strong stable and strong unstable ones as well 
as some other slow manifolds satisfying some non-resonance conditions. 

To avoid lengthening the paper, we restrict in this paper just to the one- 
dimensional manifolds (see Section [5]) , where we do not need to deal with 
resonances as it is the case in higher dimensions. We think that, considering 
this particular case, we can state in a more clear and simpler way the main 
idea behind the algorithms. We will come back to the study of higher 
dimensional manifolds in future work. 
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3.5.1. Invariant manifolds of rank 1. We once again use a parameteriza- 
tion to describe the whiskers. This amounts to finding a solution u of the 
equations of motion under the form 

u (") = W(oon,fi n s) 

in the discrete time case and 

u(t) = W(u}t,se^) 

in the continuous time case, where W : T £ x (V C M. d ~ e ) — > M. and fx € M.. 
The function W has then to satisfy the following invariance equations 

F(W(6,s)) = W{9 + u;,fis), 

d (17) 
d u W{9, s) + »Sq-W{9, s) = (Xo W){0, s), 

for the case of maps and flows, respectively. 

Note that equations (|17p imply that in variables (9, s) the motion on the 
torus consists of a rigid rotation of frequency u whereas the motion on the 
whiskers consists of a contraction (or an expansion) by a constant /x (e^ in 
the case of flows). We call contractive the situation < 1 for maps (or 
fi < for flows). We call expansive the case when > 1 for maps (or ji > 
for flows). Note that if W(9,s) satisfies CE7|) then W(9,0) is a solution of 
the invariance equations dl|) or (ITT|) . 

As in the case of invariant tori, it will be convenient to consider the 
following modified invariance equations 

F(W(9, s)) = W(9 + u, lis) + {J{K )- l DK Q ) o T^X, 

d ( 18 ) 
d 0J W(9,s)+fis—W(e,s) = (XoW)(9,s) + J(K y\DX o K )X, 

OS 

where Kq is, as before, a given approximate solution of the equations ([5]) 
and (fTT]) . respectively. 

3.5.2. Uniqueness of solutions of the invariance equation for whiskers. The 
solutions of equations (fTT]) are not unique. Indeed, if W(9, s) is a solution, 
for any a € T , b G M., we have that W(9, s) = W(0 + o~, sb) is also a solution. 
This non-uniqueness of the problem can be removed by supplementing the 
invariance equation with a normalization condition. 

Some suitable normalization conditions (in the case of maps) that make 
the solutions unique are 



W(9,0)-9 = 0, 

DF(W(6, 0))D 2 W(9, 0) = nD 2 W(9, 0), ^ 19 ^ 
||W(-,0)||=p 

where D 2 W denotes the derivative with respect to the second argument, 
p > is any arbitratrily chosen number and |.| stands for a suitable norm. 
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The fact that the solutions of ([8]) supplemented by (fTUl) are locally unique 
is proved in jFdlLS09a]. In this paper, we will see that these normalizations 
uniquely determine the Taylor expansions (in s) of the function W whenever 
the first term W\{9) = D2W(9,0) is fixed, and we will present algorithms 
to perform these computations. 

The first equation in (119j) amounts to choosing the origin of coordinates 
in the parameterization of the torus and, therefore eliminates the ambiguity 
corresponding to a. (Check how does (|19p change when we choose a). 

The second equation in (fl9l) indicates that W% (9) is chosen to be a vector 
in the hyperbolic direction. We furthermore require that we have chosen 
the coordinate so that it is an eigenvector of the expanding/contracting 
direction. 

The third equation in (|19p chooses the eigenvalue. Equivalently, it fixes 
the scale in the variables s. Observe that, setting b amounts to multiplying 
W\ by b. Hence, setting the norm of DW sets the b. 

From the mathematical point of view, all choices of p are equivalent. 
Nevertheless, from the numerical point of view, it is highly advantageous to 
choose 1 1 Wi|| so that the numerical coefficients of the expansion (in s) of W 
have norms that neither grow nor decrease fast. This makes the computation 
more immune to round off error since it becomes more important when we 
add numbers of very different sizes. 

3.6. Fourier- Taylor discretization. One of the ingredients of algorithms 
to solve the functional equations is to consider discretizations of functions 
one searches for. 

In this section, we introduce the discretizations we will use. Roughly, 
for periodic functions, we will use both a Fourier series discretization and 
a real discretization on a grid. We will show that the Newton step can be 
decomposed into substeps which require only O(N) operations in either of 
the representations. Of course, one can switch between both representa- 
tions using 0(N]n(N)) operations using FFT algorithms. For the study of 
invariant manifolds, we will use Taylor series in the real variables. 

3.6.1. Fourier series discretization. Since we are seeking functions K which 
are periodic in the angle variable 9, it is natural to discretize them retaining 
a finite number of their Fourier coefficients 

K(9)= c ^ 2mk -\ (20) 

where 

Q N = |jfe G 1} | < iv} . 

Since we will deal with real-valued functions, we have ct = c~k and one can 
just consider the cosine and sine Fourier series, 

K(9)=a + ^ a k cos(2irk-9) + b k sm(2-jrk-9). (21) 

keO N 
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These Fourier discretizations have a very long history going back to classi- 
cal astronomy, but have become much more widely used with computers and 
go under different name such as "automatic differentiation". The manipula- 
tion of these polynomials are reviewed in [Knu97j. A recent review of their 
applications in dynamics - including implementation issues and examples - 
is |Har08| . 

The main shortcoming of Fourier series discretization of a function is that 
they are not adaptative and that for discontinuous functions, they converge 
very slowly and not uniformly. These shortcomings are however not very 
serious for our applications. 

Since the tori are invariant under rigid rotations, they tend to be very 
homogeneous, so that adaptativity is not a great advantage. Also, it is 
known [FdlLS09a] that if tori are C r for sufficiently large r, they are in fact 
analytic. 

The fact that the Fourier series converge slowly for functions with dis- 
continuities is a slight problem if one wants to compute tori close to the 
breakdown of analyticity, when the tori transform into Aubry-Mather ob- 
jects. Of course, when they are far from breakdown - as it happens in many 
interesting problems in celestial mechanics - the Fourier coeffients converge 
very fast. To perform calculations close to breakdown, the a posteriori theo- 
rems in [FdlLS09a] prove invaluable help to have confidence in the computed 
objects. 



3.6.2. Fourier vs grid representation. Another representation of the function 
K is to store the values in a regularly space grid. For functions of t variables, 
we see that if we want to use N variables, we can store either the Fourier 
coefficients of index up to 0(N 1 ^) or the values on a grid of step 0(N~ 1 / e ). 

Some operations are very fast on the real space variables, for example 
multiplication of functions (it suffices to multiply values at the points of 
the grid). Also, the evaluation of F o K is very fast if we discretize on a 
grid (we just need to evaluate the function F for each of the points on the 
grid). Other operations are fast in Fourier representation. For example, it 
is fast to shift the functions, to take derivatives and, as we will see in (13. 7h . 
to solve cohomology equations. Hence, our iterative step will consist in the 
application of several operations, all of which being fast - O(N) - either in 
Fourier mode representation or in a grid representation. Of course, using the 
Fast Fourier Transform, we can pass from a grid representation to Fourier 
coefficients or viceversa in 0(N In N) operations. There are extremely ef- 
ficient implementations of the FFT algorithm that take into account not 
only operation counts but also several other characteristics (memory access, 
cache, etc.) of modern computers. 
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3.6.3. Fourier- Taylor series. For the computation of whiskers of invariant 
tori, we will use Fourier- Taylor expansions of the form 

oo 

W{e,s) = Y J W n (9)s n , (22) 

n=0 

where W n are 1-periodic functions in 6 which we will approximate using 
Fourier series (f20|) . 

To manipulate this type of series we will use the so called automatic dif- 
ferentiation algorithms (see |Knu9 7].[Har08]). For the basic algebraic opera- 
tions and the elementary transcendental functions (exp, sin, cos, log, power, 
etc.), they provide an expression for the Taylor coefficients of the result in 
terms of the coefficients of each of the terms. 

3.7. Cohomology equations and Fourier discretization. In the New- 
ton step to construct KAM tori, one faces solving cohomology equations, 
that is, given a periodic (on T^) function n, we want to find another periodic 
function ip solving (the first equation is a small divisor equation for flows 
and the second one for maps) 

cW = V, 

( 23 ) 

ip - <p o T w = n. 



As it is well known, equations (|23j) have a solution provided that 

Vo = f V = 0, (24) 

and that ui is Diophantine in the appropriate sense. The Fourier coefficients 
0k of the solution ip of (|25j) are then given respectevely by 

Vk 



2iriuj • k 

(25) 

„ _ Vk 

™ ^ _ p2-wik-u) ' 

where % are the Fourier coefficients of the function r\. 

Notice that the solution cp is unique up to the addition of a constant (the 
average (p® of ip is arbitrary). 

Equations (|23p and their solutions (|25p are very standard in KAM theory 
(see the exposition in [LlaOl]). Very detailed estimates can be found in 
[Riis75| , when oj is Diophantine (which is our case) . 

4. Fast Newton methods for (possibly) whiskered tori 

In this section we develop an efficient Newton method to solve the invari- 
ance equations (|8l)- (fTT1) and (fT0l) - (fT3l) . We mainly focus on the case of maps 
(the case for vector fields being similar is described in the appendices). 
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We emphasize that the algorithm applies both to whiskered tori and to 
Lagrangian tori. Indeed, the case of Lagrangian tori is simpler. The hyper- 
bolic part of the Lagrangian tori is just empty so that we do not need to 
compute the splittings. We refer to Algorithm [T] and Remark El 

We will assume that the motion on the torus is a rigid rotation with a 
Diophantine frequency uo € K . As we have already shown, the invariance 
implies that the vectors in the range of DK are invariant under DF. The 
preservation of the symplectic structure implies that the vectors in the range 
of (J o K)~ l DK grow at most polynomially under iteration. We note also 
that tori with an irrational rotation are co-isotropic, (DK) T (J o K)~ 1 DK = 
0, i.e. 

RangeDKnRange(JoK)- 1 DK = {0}, (26) 

and therefore dim(Range DK © Range (J o K)~ 1 DK) = 21. Hence, at any 
point of the invariant torus of dimension I with motion conjugate to a ro- 
tation, we can find a 2£-dimensional space of vectors that grow at most 
polynomially under iteration. As it is shown in [FdlLS09b], approximately 
invariant tori are approximately co-isotropic and the transversality (|26p also 
holds. 

The tori that we will consider are as hyperbolic as possible, given the 
previous argument. That is, we will assume that there exist directions that 
contract exponentially in the past or in the future, which span the comple- 
ment of the tangent to the torus and its symplectic conjugate. 

We will consider tori that have a hyperbolic splitting 

T K {6)M = £^{6) © £ k(9) © £ K{ey ( 27 ) 

such that there exist < li\,Hi < 1, Hz > 1 satisfying //1/X3 < 1, /X2M3 < 1 
and C > such that for all n > 1 and 6 G 



v 6 £|r(0) 



\M(n,9)v\ < Cni\v\ 
\M{n,6)v\ < 
\M{n,6)v\ < Cn%\v\ 



where A4(n,9) is the cocycle with generator Z(6) 
quency u, i.e. M : Z x T f ^ GL(2d,R) is given by 



Vn > 1 

Vn < 1 (28) 
Vn G Z 

= DF(K(9)) and fre- 



M(n,6) 



Z(9 + (n-l)ui)---Z(9) 
Id 



n > 1, 
n = 0, 



We will also assume that 
dim£ c K{9) = 2£, 



+ (n + l)w)---Z- L (9) n<l 



dim£ s K{e) 



d 



(29) 



(30) 



The assumption (I30p implies that the only non-hyperbolic directions are 
those spanned by the tangent to the torus and its symplectic conjugate, 
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that is, there are no elliptic directions except those that are forced by the 
symplectic structure and the fact that the motion on the torus is a rotation. 

We associate to the splitting (|27|) the projections IL c K ,gy ^-KfO) anc ^ ^k(S) 
over the invariant spaces £^(m, ^k(Q) anc ^ ^K(8Y 

It is important to note that since F is symplectic (i.e. F*Q = O), for all 
n > 1 and n < — 1 

tt(u,v) = n(DF n u,DF n v), 
so that, if u, v have rates of decrease, by taking limits in the appropriate 
direction we obtain that O is zero. That is, we get 

n(s s ,s s ) = o, n{£ u ,s u ) = o, 
n(£ c ,s s ) = o, n(£ c ,£ u ) = 0. 

Therefore, we have 



(J(K(8))Y 


-\ F c 


- C K(8) 


(j(K(9))y 


-l F s 


cu 

- C K(6) 


(j(K(e))y 


-1 CU 

C K(8) 


— C K(9) 



In [HdlLS09], we provide a method to compute the rank-1 bundles by 
iterating the cocycle. Of course, once we have computed the vector spanning 
the rank-1 (un)stable bundle it is very easy to obtain the projections. In 
Section [4.31 we discuss an alternative to compute the projections by means 
of a Newton method. In this case we do not need to assume that the bundle 
is 1-dimensional. 

4.1. General strategy of the Newton method to solve the invari- 
ance equation. In this section we will design a Newton method to solve 
the invariance equation ([8]) and the modified one (flOl) . and discuss several 
algorithms to deal with the linearized equations. 

We first define the following concept of approximate solution. 

Definition 1. We say that K (resp. (K, A)) is an approximate solution of 
equation ((8]) (resp. CE])) if 

FoK-KoT ul = E, 

(31) 

(resp. FoK-KoT u - ((J o K^DKq) o T u \ = E), 
where E is small. 

The Newton method consists in computing A in such a way that setting 
K K + A and expanding the LHS of (i3T|) in A up to order ||A|| 2 , it 
cancels the error term E. 

Remark 5. Throughout the paper, we are going to denote ||.|| some norms 
in functional spaces without specifying what they are exactly. We refer the 
reader to [LGJV05, FdlLS09a], where the whole theory is developed and the 
convergence of the algorithms is proved. Recall that one of the key ideas of 
KAM theory is that the norms are modified at each step. 
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Performing a straightforward calculation, we obtain that the Newton pro- 
cedure to solve equation (jS]) and (fTT|) . given an approximate solution K, 
consists in finding A satisfying 

(DF o K)A - A o T w = -E. (32) 



For the modified invariance equation (|10p . given an approximate solution 
(K, A), the Newton method consists in looking for (A, S) in such a way that 
K + A and A + 5 eliminate the error in first order. The linearized equation 
in this case is 

(DF o K)A - AoT,-((Jo Ko^DKo) o TJ = -E, (33) 

where one can take Kq = K. 

As it is well known, the Newton method converges quadratically in ||_E7|| 
and the error E at step K + A is such that 

\\E\\ < C\\E\\ 2 

where E is the error at the previous step. 

In order to solve the linearized equations (I32D and (I33D . we will first 
project them on the invariant subspaces £ c ,£ u and 8 s , and then solve an 
equation for each subspace. 

Thus, let us denote 

A^(0) = n^;A(#), 

(34) 

E^{6) = U s ^ +w) E(9), 

such that A(6>) = A s (6>) + A c (6>) + A"(6»). Then, by the invariant properties 
of the splitting, the linearized equations for the Newton method (|32l) and 
(|33l) split into: 

DF(K(6))A C {0) - A c o T u (0) = -E c {9), 
DF(K(6))A S {9) - A s o T u {6) = -E s (8), (35) 
DF{K{9))A U {0) -A u o T u {6) = -E u (6), 

and 

DF{K(9))A C (9) - A c o T u {0) + U c K{0+uj) (J o K (9 + u^DK^ + u)8 = -E c (8), 

DF{K(9))A S (9) - A s o T u {0) + U s K{s+Lu) (J o K o (0 + u))- x DK Q {8 + oj)5 = -E s (8), 

DF{K(9))A U (9) - A u o T w (0) + n^ (e+aj) ( J o K o (0 + u))- l DK Q {8 + u)6 = -E u {8). 

(36) 

Notice that once 5 is obtained, the equations (j36|) on the hyperbolic spaces 
reduce to equations of the form ([35]) . More precisely, 

DF(K(8))A S > U (8) - A S ' M o T u (0) = -E s > u {9) (37) 
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where 

= E s,u (e) + n^ +£j) (J o K (9 + u^DKoie + u)5. 

Equations ([35]) and ([36]) for the stable and unstable spaces can be solved 
iteratively using the contraction properties of the cocycles on the hyperbolic 
spaces given in (j28]h Indeed, a solution for equations ([37|) is given explicitly 
by 

oo 

A s (9) = E s oT^ w (9)+Y,{DFoKoT^{9)x- ■ ■xDFoKoT- ku (6))(E'oT Hk+1)l 

k=l 

(38) 

for the stable equation, and 

oo 

= - ^(DF' 1 K $) x ' ' ' x DF ~ X K ° T U#))(£ U ° T ^(^)) ( 39 ) 

fc=0 

for the unstable direction. Of course, the contraction of the cocycles guar- 
antees the uniform convergence of these series. 

The algorithms presented in Appendix [A] allow us to compute the solu- 
tions A S ' M of equations (|3T[) efficiently. 

In Section 14.21 we discuss how to solve equations (|35p and (|36p in the 
center direction. 

Hence, the Newton step of the algorithm for whiskered tori that we sum- 
marize here will be obtained by combining several algorithms. 

Algorithm 1. Consider given F, u, Kq and an approximate solution K 
(resp. K,X), perform the following operations: 

A) Compute the invariant splittings and the projections associated to 
the cocycle Z(9) = DF o K{9) and uj using the algorithms described 
in Section\J^\ (or in |HdlLS09] ). 

B) Project the linearized equation to the hyperbolic space and use the 
algorithms described in Appendix\A\ to obtain A s ' u . 

C) Project the linearized equation on the center subspace and use the 
Algorithm^ in Section \4.2\ to obtain A c and S. 

D) Set K + A s + A u + A c -»• K and A + 5 -> A 

Of course, since this is a Newton step, it will have to be iterated repeatedly 
until one reaches solutions up to a small tolerance error. 

We will start by some remarks on the different steps of Algorithm [1] and, 
later, we will provide many more details on them. 

Remark 6. It is important to remark that the above Algorithm [T] also applies 
to the case of Lagrangian tori. It suffices to remark that in that, case, the 
center space is the whole manifold, so that there is no need to compute the 
splitting. Hence, for Lagrangian tori, the steps A) and B) of Algorithm Q] 
are trivial and do not need any work. 
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Remark 7. The main issue of the Newton method is that it needs a good ini- 
tial guess to start the iteration. Any reasonable algorithm can be used as an 
input to the Newton method. Indeed, our problems have enough structure 
so that one can use Lindstedt series, variational methods, approximation 
by periodic orbits, frequency methods, besides the customary continuation 
methods. 

Remark 8. As we have mentioned in Remark [3j the solutions of ([8]) and (llip 
are not unique. Therefore, in order to avoid dealing with non-invertible ma- 
trices in the Newton procedure, we will impose the normalization condition 

/ (K(e)-I(6))-u l = 

where {ui}f =1 is a basis for Range(I) (L being the dimension) and I is a 
linear function introduced in (|15p . 

4.2. Fast Newton method for (whiskered) tori: the center direc- 
tions. We present here the Newton method to solve the equations on the 
center subspace in the case of maps. 

For Lagrangian tori, the hyperbolic directions are empty and the study 
of the center direction is the only component which is needed. Hence, the 
algorithms discussed in this section allow to solve, in particular, equations 
(|32p and (|33p in the case of Lagrangian tori. For a discussion of the center 
equations for Hamiltonian flows, we refer the reader to Appendix IB"1 

The key observation is that the linearized Newton equations (|32p and 
(f33"P are closely related to the dynamics and therefore, we can use geomet- 
ric identities to find a linear change of variables that reduces the Newton 
equations to upper diagonal difference equations with constant coefficients. 
This phenomenon is often called "automatic reducibility" . 

The idea is stated in the following proposition: 

Proposition 2 (Automatic reducibility, sec FdlLS09b~l IFdlLS09a] ). Given 
an approximation K of the invariance equation as in (|3ip . denote 



a{9) 


= DK(6) 


N{9) 


= ([a{6)]M0))~ l 


m 


= a{6)N(9) 


7(0) 


= {joK{6))- i m 



(40) 



and form the following matrix 

M(9) = [a(6) | 7(^)1, (41) 

where by [• \ ■] we denote the 2d x U matrix obtained by juxtaposing the two 
2d x £ matrices that are in the arguments. 
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Then, we have 

(DF o K(9))M(9) = M(9 + u) ^ A W) + E(9) (42) 

where 

A{9) = (3(9 + ^[(DF o K(9)) 1 {9) - 7 (0 + u)], (43) 
and \\E\\ < \\DE\\ in the case of (|32p or \\E\\ < \\DE\\ + |A| in the case of 

Remark 9. If the symplectic structure is almost-complex (i.e. J 2 = — Id), 
we have that 

(3(9 + ^^(9 + ^) = 0, 
since the torus is isotropic. Then A(9) has a simpler expression given by 

A{9) = 13(9 + uj) l (DF o K)(9)j(9). 

Once again, we omit the definition of the norms used in the bounds for E. 
For these precisions, we refer to the paper [FdlLS09a], where the convergence 
of the algorithm is established. 

It is interesting to pay attention to the geometric interpretation of the 
identity (I42p . Note that, taking derivatives with respect to 9 in (|3ip . we 
obtain that 

(DF o K)DK - DK oT u = DE, 

which means that the vectors DK are invariant under DFoK (up to a certain 
error). Moreover, (J o K)~ 1 DKN are the symplectic conjugate vectors of 
DK, so that the preservation of the symplectic form clearly implies (I42D . 
The geometric interpretation of the matrix A(9) is a shear flow near the 
approximately invariant torus. See Figured! 

DF(K(9))v(9) 




To be able to use the change of unknowns via the matrix M previously 
introduced on the center subspace, one has to ensure that one can identify 
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the center space £j((d) w ith the range of M. This is proved in |FdlLS09aj to 
which we refer. 

For our purposes it is important to compute not just the invariant spaces, 
but also the projections over invariant subspaces. Knowing one invariant 
subspace is not enough to compute the projection, since it also depends on 
the complementary space chosen. 

In the following, we will see that the result stated in Proposition [2] allows 
us to design a very efficient algorithm for the Newton step. 

Notice first that if we change the unknowns A = MW in (|32p and (|33p 
and we use (1421) we obtain 



M(6 + u) A jP) W{9) - M(9 + oj)W{6 + u) 

V u id / (44) 

- (J(K (8 + u)))- l DK (6 + u)S = -E{9) 



Of course, the term involving S has to be omitted when considering (|32p . 

Multiplying flM]) by M(9 + uj) 1 - J{K(9 + u)) and using the invertibility of 
the matrix M(0+uj) l J(K{0 + oj))M{0 + uj) (see |FdlLS09bllFdlLS09a] \ wc 
are left with the system of equations 

Wi(0) + A{9)W 2 {9) - B x {e)8 - Wi{9 + oj) = -E^O) 

(45) 

W 2 {6) -W 2 (6 + u)- B 2 (9)5 = -E 2 (9) 

where 

E{8) = (M(9 + uj) l J{K{9 + lo))M(9 + lo))- 1 M(9 + co) ± J(K(9 + u))E{0) 

B{9) = {(M ± J(K)M)- 1 M ± J(K)(J(K ))- l DK }oT UJ (9) 

and the subindices i = l,2 indicate symplectic coordinates. 

When K is close to Kq, we expect that B 2 is close to the ^-dimensional 
identity matrix and B\ is small. 

The next step is to solve equations (f45l) for W (and 5). Equations (l45l 
are equations of the form considered in (I23p and they can be solved very 
efficiently in Fourier space. 

More precisely, the second equation of (j45[) is uncoupled from the first 
one and allows us to determine W 2 (up to a constant) and 5. The role of the 
parameter 5 is now clear. It allows us to adjust some global averages that 
we need to be able to solve equations ([45]) , Indeed, we choose 5 so that the 
term B 2 {9)5 — E 2 has zero average (which is a necessary condition to solve 
small divisor equations as described in Section [3 .7p . This allows us to solve 
equation ([23]) for W 2 . We then denote 

W 2 (9) = W 2 (9) + W 2 

where W 2 (9) has average zero and W 2 S E. 
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Once we have W2, we can substitute W 2 in the first equation. We get W 2 
imposing that the average of 

B x (9)6 - A{9)W 2 {9) - A{9)W 2 {9) - E^O) 

is zero and then we can find W\ up to a constant according to (f25|) . 

We therefore have the following algorithm to solve ([3]) in the center di- 
rection, 

Algorithm 3 ( Newton step in the center direction). Consider given F, 
uj, Kq and an approximate solution K (resp. K,\). Perform the following 
calculations 

1. (1.1) Compute F oK 

(1.2) Compute K o T w 

(1.3) Compute the invariant projections, rP,n",II c . 

2. Set E c = U C (F o K - K o T u ) (resp. set E c = Il c (F oK-KoT u - 
(JoKo^DKoX)) 

3. Following (@EJ) 

(3.1) Compute a(6) = DK(0) 

(3.2) Compute N{9) = ([a(9)] ± a(9)) _1 

(3.3) Compute /3(9) = a(9)N(9) 

(3.4) Compute j(9) = (J(K(9)))- 1 /3(9) 

(3.5) Compute M{9) = [a(9) | 7 (0)j 

(3.6) Compute M(9 + uj) 

(3.7) Compute (M(9 + u) ± J(K(9 + w))M(0 + w))" 1 

(3.8) Compute E{9) = {M{9 + u) x J(K(9 + uj))M{9 + w))- 1 .^). 
We denote E±,E 2 the components of E along DK and along 
J- l DK. 

(3.9) Compute 

A{9) = f3{9 + uj) l [{DF o K(9)) 1 (9) - 7 (0 + u)] 

as indicated in (|43p 

4. (4.1) Solve for W2 satisfying 



W 2 -W 2 oT u) = -E 2 - f E 2 

Jt 1 



( resp. 

(4.1') Solve for 5 such that 







[ Ez- 


[ B 2 


Jr* 





■ 

(4.2') Solve for W 2 satisfying 

W 2 -W 2 oT UJ = -E 2 + B 2 5 

Set W 2 such that the average is 0.) 
5. (5.1) Compute A{9)W 2 {9) 
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A{9) 



Wo 



A{9) 



W-2 



(5.2) Solve for W 2 satisfying 
0=1 E 1 (9)+ [ A{9)W 2 {9) + 

(5.3) Find W\ solving 

Wi - Wi oT u = -Ex - A(W 2 +W 2 

Normalize it so that jj e W\ = 
( resp. 

(5.1') Compute A{9)W 2 {9) 
(5.2') Solve for W 2 satisfying 

0=1 E x {9)- I B 1 (9)5 + I A(9)W 2 {9) + 

Jje JT l JT e 

(5.3') Find W\ solving 

W 1 -W 1 oT u = -E x - A(W 2 + W 2 ) + B x 5 

Normalize it so that jj e W\ = 0.) 

6. The improved K is K{9) + M(9)W{9) 
(resp. the improved A is A + 5). 

Notice that steps (1.2), (3.1), (3.6), (4.1) (resp. (4.2')), (5.3) (resp. (5.3')) 
in Algorithm [3] are diagonal in Fourier series, whereas the other steps are 
diagonal in the real space representation. Note also that the algorithm only 
stores vectors whose size is of order N . 

Remark 10. Using the symplectic properties of the matrix M, step (3.8) can 
be sped up. 

When the torus is exactly invariant we have that the invariant torus is 
co-isotropic. Hence DK^JoKDK = 0. Hence, when the torus is invariant, 
we have 

M{9 + lo) ± J(K{9 + u))M{9 + u)= (^_° N ^ 

so that the inverse is easy to calculate. 

For the purposes of a Newton Method, we can use the same expression 
for the inverse in step (3.8) and still obtain a quadratically convergent algo- 
rithm. 

4.3. A Newton method to compute the projections over invariant 
subspaces. In this section we will discuss a Newton method to compute 
the projections Il c K ^y ^-jcfB) an d ^k($) associated to the linear spaces £j(/g\, 
£ S K ^ and £^(0) where K is an (approximate) invariant torus. More precisely, 
we will design a Newton method to compute ITw^ and n™,^ = H° K (g) + 
~^k(oy Similar arguments allow us to design a Newton method to compute 

U K(9) and U K(e) = W K(9) + U K(e)- Then > ° f course > n K(9) is S iven b y 

ttc TJ CS TJ CU TJ CU TJ CS 

ll K(9) — il K(e) iL K(9) — ll K(0) ll K(B) ' 
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Let us discuss first a Newton method to compute R s K rg\ and ITwm. To 
simplify notation, from now on, we will omit the dependence in K(6). 

Given a cocycle Z(9) (which in our case will be Z{9) = DF(K(Q))), we 
will look for maps IP : -> GL(2d,R) and IP" : -> GL(2d,R) satisfying 
the following equations: 



u cu (e + uj)z{e)u s (e) = o, 


(46) 


u s (e + u)z(e)u cu (9) = o, 


(47) 


u s (6) + u cu (e) = id, 


(48) 


[u s (0)} 2 = n s {6), 


(49) 


[U cu {9)] 2 = U cu {9), 


(50) 


U s (0)U cu (0) = 0, 


(51) 


U CU (9)U S (9) = 0. 


(52) 



Notice that the system of equations ([4^ - ([52]) is redundant. It is easy to 
see that equations (jSTfl) . (jSTj) and (j52j) follow from equations (f48|) and (|49|) . 
Therefore, the system of equations that needs to be solved is reduced to 
equations (l46j) - (li9j) . 

We are going to design a Newton method to solve equations (j4"6l) - ([4"7|) 
and use equations (I48p ~ (l49p as constraints. In this context, by approximate 
solution of equations (f4U|) - (}4"T|) . we mean a solution (II s , H cu ) such that 



U cu {9 + u)Z(0)U s {9) = E cu (9), (53) 

U s (9 + oj)Z(9)U cu (9) = E s (9), (54) 

U s (9) + U cu (9) =Id, (55) 

[U s (d)] 2 = U S (9). (56) 



where E % denotes the error in a certain component. Notice that the error 
in equation (|53p has components only on the center and unstable "approx- 
imated" subspaces and we denote it by E cu . The same happens with the 
equation ([54"}) but on the "approximated" stable subspace. We assume that 
E cu and E s are both small. 

As standard in the Newton method, we will look for increments A s and 
A cu in such a way that setting IP <- IP + A 3 and IP" <- IP" + A cu , 
the new projections solve equations (|46p and (|47|) up to order ||-E|| 2 where 
\\E\\ = \\E S \\ + ||-E C "|| for some norm ||.||. 

The functions A s and A cu solve the following equations 

A cu {9 + u)Z(9)U s (9) + n CM (0 + u)Z(6)A s (6) = -E cu (6) 

a s (9 + co)z(e)u cu (e) + n s (^ + u)z{9)A cu {e) = -e s {9) [ ' ' 1 

with the constraints 



A s (9) + A cu {9) = 

U S (9)A S (9) + A s (0)n s (^) = A s (9) 



(58) 
(59) 
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By equation (I58|) we only need to compute A s since A cu = —A s . We now 
work out equations (|57l) . (f58l) and (i59j) so that we can find A s . 
Denote 

A s = U CU A S ^ 

so that 

A* = A* + A*„. (61) 

Then equation (I59D reads 

A*(0) + A'(0)H'(0) = A s s (0) + A^(0), (62) 

or equivalently, 

A s (e)U s (0) = A s cu (6) . (63) 

By d5SI), dMI) and dM]) we have that 

A S (9)U CU (9) = A s (9) - A S (9)U S (9) = A s (9) - A s cu (6) = A s s (6). (64) 

Now, using (f58j) . equations (f57|) transform to 

- A s (fl + oj)Z(9)U s (9) + n cu (# + w)Z(#)A s (#) = -E cu (9), 
A s (9 + u)Z(9)U cu {9)-U s {9 + uj)Z{9)A s {9) = -E s {9). l05) 

Denoting 

at s (0) = n s (e + w)z(0)n s (0), 
at cu (0) = n c "(0 + w )z(#)n cu (#), 

and using that n cu (<9 + oj)Z{9)W(9) and n s (6> + u)Z(9)U cu (9) are small by 
([55 ]) -([5I |l and U s (9) + n c "(6>) = Id by ([55]). it is enough for the Newton 
method to solve for A s satisfying the following equations 

- A s (9 + u)U s (9 + u)N 3 (6) + N cu (9)n cu (9)A s (9) = -E cu (9), 

A s {9 + u)U cu (9 + u)N cu {9) - N S (9)U S {9)A S (9) = -E s {9). l0<>) 

Finally, by expressions (|63p and ()64|) and taking into account the notations 
introduced in ([60]) . equations ([66]) read 



-A s cu (9 + a,)jV.(0) + N^A^O) = -E cu (9), (67) 
A*(0 + u)N cu (9) - N S (9)A S S (9) = -E s (9). (68) 

In Appendix^] we discussed how to solve efficiently equations of the form 
dSTD-dSH]). Notice that they are of the form {97]) for A(9) = N cu (9), B(9) = 
N s {9) and r/(6>) = -E cu (9) in the case of equation {67| and A(9) = N s (6), 
B{9) = N cu {9) and r](9) = +E S (9) in the case of equation (jM]). Further- 
more, ||iV s || < 1 and H-A^H < 1. Hence, they can be solved iteratively using 
the fast iterative algorithms described in Appendix lAl 
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The explicit expressions for Af, u and A* are 



cu 



A s cu (9) = - N- U 1 (9)E CU (9) + N- u \6) x • - - x 

n=l (69) 

N^{6 + nu)E cu {9 + nu)N s (9 + (n - l)w) X •• • X N s (6) 

and 

oo 

A s s (9) = E s (9 - u)N- u l {9 - W ) + £ JV.(0 - u) x • • • x 

n=l 

N s (9 - (n + l)u;)£ s (0 - (n + l)u)N£(6 - (n + l)w) x • • • x AT" 1 (5 - w). 

(70) 

Remark 11. Notice that by the way N cu (9) is defined, it is a matrix which 
does not have full rank. Therefore, we understand N^(9) as the "pseudoin- 
verse" matrix. 

Finally, let us check that A s = A^ u + also satisfies the constraints. In 
order to check that constraint ([59]) . which is equivalent to ([63]) . is satisfied 
we will use the expressions (I69p and flTOH . Notice first that 



N S (9)U S (9) = N S (9) (71) 

and 

N- u 1 (9-u)U s (9) = 0. (72) 

Moreover, from (153D and using (156p one can see that 

E CU (9)U S (9) = U cu (9 + oj)Z(9)[U s (9)} 2 = E cu {9) . (73) 

Then, from expressions ([69]) and (fT0|) and the above expression ([71]) . ([72]) 
and ([73]) . it is clear that 

A S (9)U S (9) = A S S (9)U S (9) + A S CU (9)U S (0) = + A* M , 



hence, constraint ([63]) is satisfied. 

Now, using equation ((SSJ we get A c "(6») = -(A*(0) + A^(6»)) and the 
improved projections are 

n s (9) = n s (9) + A s s (9) + A s cu (9) 
U cu (9) = U cu (9) + A cu {9). 

The new error for equations (I46p and (I47D is now < Clji?!! 2 where 
| (E 1 1 1 = \\E CU \\ + \\E S \\. Of course equation ((3HD is clearly satisfied but ([4"9]) 
is satisfied up to an error which is quadratic in \\E\\. However it is easy to 
get an exact solution for (I49p and the correction is quadratic in A s (and 
therefore in A cu ). To do so, we just take the SVD decomposition of II s and 
we set the values in the singular value decomposition to be either 1 or 0. 
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In this way we obtain new projections II new and n^" w = Id — II new satis- 
fying 

HIT - n s ii < ||A S || 2 

1 1 -"new II ^ II I 

MTTCti ttCU | . II A CU ||2 

Harlow — II ^ 11^ II > 

so that the error for equations (|46p and (|47p is still quadratic in \\E\\. More- 
over, they satisfy equations (|4"9j) and, of course, (jUJ) exactly. 

Hence, setting II s <— II new and II C " <— n^ w we can repeat the procedure 
described in this section and perform another Newton step. 

Consequently, the algorithm of the Newton method to compute the pro- 
jections is: 

Algorithm 4 (Computation of the projections by a Newton method). Con- 
sider given F,K,u and an approximate solution (II s , IF") of equations (|46[) - 
(|47|) . Perform the following calculations: 

1. Compute Z{8) = DF o K{9) 

2. (2.1) Compute E cu {8) = U cu (9 + u)Z(9)W{9) 
(2.2) Compute E s {6) = W{8 + uj)Z(9)W u (9) 

3. (3.1) Compute N s (8) = U s (8 + u)Z{8)W(9) 
(3.2) Compute N cu (9) = U cu (8 + oj) Z(8)U CU (9) 

4. (4.1) Solve for A s s satisfying 

N S (0)A S S (9) ~ A*(0 + co)N cu (0) = E s {9) 
(4.2) Solve for A s cu satisfying 

N CU (8)A S CU (8) - A cu (6 + u)N.(P) = -E cu (8) 

5. (5.1) Compute U s (9) = U s (9) + A s s {8) + A cu (9). _ 

(5.2) Compute the SVD decomposition ofW{9): IF(0) = U{9)^{9)V L {9). 

(5.3) Set the values in X(#) equal to the closer integer (which will be 
either or 1). 

(5.4) Recompute W(6) = U{9)^(9)V ± {9). 

6. Set IT -> IT 

id -n s n cM 

and iterate the procedure. 

Notice that the matrix multiplication is diagonal in real space represen- 
tation, whereas the phase shift is diagonal in Fourier space. A discussion on 
how to perform step 4 efficiently is given in Appendix IA1 

5. Computation of rank-1 whiskers of an invariant torus 

In this section, we present algorithms to compute the whiskers associated 
to an invariant torus, that is the invariant manifolds that contain the torus 
and are tangent to the invariant bundles. 

For the sake of simplicity and in order to state in a clear way the main 
idea behind the methods we will only discuss the case when the invariant 
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whiskers are one-dimensional (i.e. d— £ = 1). The same idea can be extended 
to compute invariant manifolds of any rank. However, there are several new 
phenomena (resonances) that can appear and need to be discussed. We plan 
to come back to this issue in the future. 

As we already mentioned in Section 13.5.11 we will look for the whiskers 
by finding a parameterization for them, so we will search for a function 
W :T e x (V C R d ~ e ) -> M and a scalar n satisfying equation (fT7|) . 

We will consider three different methods to solve equation (|17|) . We will 
first discuss the order by order method. The other two methods are based 
on the philosophy of quasi-Newton methods. Using the phenomenon of 
"automatic reducibility" , we are able to design an efficient Newton method. 
The first method allows to compute simultaneously the invariant tori and 
the whiskers, whereas the second one assumes that the invariant tori and 
the tangent bundles are already known. 

We detail only the case of maps because the same ideas work for the case 
of vector fields and refer the reader to the appendices for the case of flows. 

Similar algorithms were developed and implemented in |HL06bl [HL07] for 
the slightly simpler case of quasi-periodic maps. 

5.1. The order by order method. In this section we adapt the parame- 
terization method introduced in [CFL05J. The convergence of the Fourier- 
Taylor series in this paper can be easily adapted to the present case. We 
focus on the case of maps and refer the reader to Appendix [O for the case 
of flows. 

We will find a solution (W, (J,) of the invariance equation (|17j) discretizing 
it in Fourier- Taylor series. Hence, we will look for W as a power series 

oo 

W(e,s) = J2Wn(0)s n , (74) 

n=0 

and match similar coefficients in s n on both sides of equation (|17p . 
For n = 0, we obtain 

F(Wo(9)) = W (9 + lo), (75) 

which is equation ([8]) for the invariant torus. Therefore, we have Wq{Q) = 
K(6), where K is a parametrization of the invariant torus. 
For n = 1, we obtain 

DFoK(6)W 1 (e) = W 1 {6 + Lu)v, (76) 

so that W\{6) is an eigenfunction with eigenvalue \x of the operator M(l, 9) 
defined in equation (|29l) . 

Equation (|76l) states that the bundle spanned by W\ is invariant for 
the linearization of F, and the dynamics on it is reduced to a contrac- 
tion/expansion by a constant [i. Therefore, one can compute W\ and \i 
using the algorithms given in Section 14.31 
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Remark 12. Notice that if W\{9) is a solution of equation (j76"j) . then bWi(9), 
for any b € M, is also a solution. See Section 13.5.21 for a discussion on how 
to choose b. 

For n > 2, we obtain 

L>F o = W" n (0 + + i? n [W , . . . , W„_i](0), (77) 

where i? n is an explicit polynomial in Wo, • • • , W n -\ whose coefficients are 
derivatives of F evaluated at Wo. 

Equation (1771) can be solved provided that fi is such that fj, n is not in 
the spectrum of the operator A4(l,9). This condition is clearly satisfied in 
the case of (un)stable bundles which are one-dimensional but it can also be 
satisfied by other bundles. 

Equation (177p can be solved using the large matrix method. It consists on 
considering a discretization of equation (|77l) using Fourier series and reduce 
the problem to solving a linear system in Fourier space, where the unknowns 
are the Fourier coefficients of the matrix W n . 

There are also efficient algorithms which are variants of the methods de- 
voted in the previous sections. The equation (j77|) is equivalent to 

W n (9) = (DF o K(9))~ 1 [v n W n (6 + co)R n [W , . . . , W n _i](0)] , 

which, for large enough n is a contraction, so that we can apply the fast 
methods of Section IA.11 In particular Algorithm IA.81 In the case that the 
stable and unstable directions are one dimensional - which is the one we 
discuss in this paper - this is enough (remember that we always have n > 2. 
When the bundles are higher dimensional, we may need to find a splitting 
corresponding to the cocycle generated by Z(6) = (DF o K(9))~ 1 fj, n . 

Remark 13. Notice that once Wq(9) and W\(9) are fixed, the solution W n {9) 
for n > 2 of equation (J77J is uniquely determined. It is then clear that any 
analytic solution is unique. The existence of analytic solutions is discussed 
in [CFL05j . 

Remark 14. Notice that the equations to compute the new term W n do not 
involve small divisors. 

Remark 15. In this case, we have not considered the modified equation (I18p 
with the counterterm, because for this method there are no obstructions to 
deal with as it is in the case of the Newton method, see Remark 1161 Indeed, 
the vanishing Lemma in [FdlLS09a] guarantees that for exact symplectic 
maps A = in (|18p once the Newton method has converged. 

5.2. A Newton method to compute simultaneously the invariant 
torus and the whiskers. In this Section we present an algorithm to solve 
equation (|17p using a Newton method, instead of solving it step by step as 
we discussed in the previous section. As before, we only deal with the case 
of maps and refer the reader to Appendix [D] for the case of flows. We do not 
prove the convergence of the algorithm here for sake of length (and purpose 
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of the paper) but it can be done using techniques similar to those developed 
in [FdlLS09bl lFdlLS09aj . 

We start with an initial approximation (W,(Jt) (resp. (W,fi, A)) satisfying 

F(W{9, s)) - W(9 + u, fis) = E{9, s) 

F(W(0, s)) - W(9 + u, fis) - J(K (9 + uj)- 1 DK (9 + co)X = E(9, s) ^ 
and we look for an improved solution 

W <- W + A 
A <- X + 5 
H <- n + 5^ 
by solving the linearized equation 

[DF o W(9, s)]A(9, s) - A(9 + u, fis) - sd s W(9 + w, ns)6 

- ((J o K^DKo) o T u (0)6a = -E(9, s). ' ' 9) 



Remark 16. As in the Newton method for invariant tori, the role of the 
parameter A is to adjust some averages to solve the equations for the case 
n = 1. More precisely, 5 will be chosen in equation (|93p so that equation 
(|9ip can be solved without any obstruction. 



We will try to solve equation (|79p by discretizing it in Fourier- Taylor 
series. Notice that equation (|79p is not diagonal when discretized in Fourier- 
Taylor series because of the term DF o W. However, there is a way to make 
it diagonal using the geometric identities which are a direct generalization 
of those used for the automatic reducibility. 

We first give the idea of the automatic reducibility when W is such that 

(F o W) {9,s) = W{6 + oj, fis) . (80) 

Taking derivatives with respect to 9 and s, we see that 

DF o W(9, s)D e W(9, s) = D e W{9 + u, fis), 

DF o W(9, s)d s W{9, s) = ^d s W(9 + u, fis). 

From the above equations, we read that the quantity D$W(9,s) remains 
invariant under DF oW(9, s), whereas the vector d s W(9, s) is multiplied by 
a factor \x. 

The vectors ( J o W)" l D e WN and (J oW)~ 1 d s WN , where N and N are 
normalization matrices (see (|8ip ). are the symplectic conjugate vectors of 
DgW and d s W, respectively. The preservation of the symplectic structure 
implies that 

(DF o W(9, s))(J(W(9, s)))- l D g W(9, s)N(9, s) = 

(J(W(9 + co, fis))Y l D e W{9 + oj, fis)N(9 + to, fis) + A(9, s)D e W(9 + co, fis), 
(DF o W(9, s))(J(W(9, s)))- 1 d s W(9, s)N(9, s) = 

- (J(W(9 + uj, fis)))~ 1 d s W(9 + oj, fis)N(9 + oo, fis) + B(9, s)d s W(9 + oo, fis). 
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where A(9,s) and B(6,s) are some matrices, which will be computed as 
before. See Proposition 03 

Therefore, one can see that in the basis DgW,(J ° W)~ 1 DqWN , d s W , 
(J o W)~~ 1 d s WN, the matrix DF o W is upper triangular with constant 
coefficients on the diagonal. 

Following the proofs in [FdlLS09a] for instance, one can prove the follow- 
ing proposition, which generalizes Proposition [2] 

Proposition 5. Denote 



a(9, s) 


= D e W(9,s) 


(3(9, s) 


= d s W{9,s) 


P(9,s) 


= a(9,s)N(9,s) 


Q(9,s) 


= (3(9,s)N(9,s) 


N[0, a) 


= {a(9,s) ± a{9,s))- 1 


N(9,s) 


= ((3(9,s) ± (3(9, s))- 1 


1(9, s) 


= (J oW(9,s))- l P(9,s) 


V(9,s) 


= (J oW(9,s))- l Q(9,s) 



(81) 



and form the following matrix 

M{9,s) = [a(9,s)\ 1 (9,s)\P(9,s)\i 1 (9,s)\ 



(82) 



where we denote by [■ \ ■ \ ■ | •] the 2d x 2d matrix obtained by juxtaposing 
the two 2d x £ matrices and the two 2d x (d — I) matrices that are in the 
arguments. 
Then 



where 



(DF o W)(9, s)M(9, s) = M(9 + u, ns)R(6, s) + 0(E), 
( _. ... . \ 

FL(9,s) = 



with 



Id A(9,s) 


o 


Id 


o 


fi B(9,s) 


1/// 



(83) 



A(9, s) = P(9, s)^ [(DF o W)(9, s) 7 (0, s) - 7 (0 + to, (is)], 



B(9, s) = Q(9, s) L [(DF o W)(9, s)n(9, s) 
and E is the error in (1781) . 



J" 



-r](9 + iv,fis)}, 
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Remark 17. If the symplectic structure induces an almost-complex one (i.e. 
J 2 = — Id), we have that 

A(9, s) = P(9, s)- L (DF o W)(9, s) 7 (9, s), 

B{9, s) = Q{9, s) ± (DF o W){9, a)rj{9, s). 

Now if we change the unknowns A = MV in (|79|) and multiply by 
M (0 + u,fj,s) the LHS, by Proposition [5j we are left with the following 
system of equations 

R{9, s)V{9, s) - V(9 + u, fts) - C{9, s)5 u = -E(9, s) + sH(9, s)6, (84) 

where R(9, s) is given in ([83]) and 

C(9,s) = M' 1 (9 + to, fxs){J{K {9 + lo)))- 1 DK (9 + u), 

E(9, s) = M~ 1 (9 + u, tis)E(0, s), 

H(9,s) = M- 1 (9 + uo^s)d s W(9 + co,fis). 

We expand the terms in (|84p in power series up to some order L and 
match coefficients of the same order on both sides of the equation. We 
use subindices to denote coordinates and superindices to denote the order. 
Hence, for order s° we have 

Vf (0) - V x °(9 + uj) + A°(9)V 2 °(9) - 01(9)5^ = -S&B), (85) 

V 2 °(9) - V 2 °(9 + w ) - Cl(9)5^ = -E° 2 (9), (86) 

»V 3 °(9) - V 3 °(9 + oo) + B°(9)V 4 °(9) - 01(9)8, = -E° 3 (9), (87) 

-V2{&) ~ V£{9 + ui)- C° 4 (9)5, = -El(9). (88) 
A* 

Notice that (185p and (|86h can be solved using Algorithm [3l Hence, we 
determine V®, V 2 and <5^. Once we know <5„, we can solve uniquely for 
and V 4 in equations (|8T|) and (|88p . These equations do not have any small 
divisors nor obstructions since 7^ 1. 

For order s 1 we have 

V{{9) - nV^{9 + 00) + A\9)V 2 \9) + A\0)V 2 °(9), (89) 

= -E\{9) + 8H° 1 {9) + 5,C\{9) 
V 2 l (9) - nY}(0 + oj) = -E\(9) + 5H$(9) + 6^(9), (90) 
liVi{9) - nV£(0 + oj) + B°(9)Vi(9) + B l (9)V A °(9) (91) 

= -El(9)+8Hi(9) + SpCltf), 

-Vl(9) - fjVl(e + oj) = -E\{9) + 5Hl(9) + 5,C\{9). (92) 
A 1 



Notice that once we choose 5, equations (I89p and (I90p are uniquely solv- 
able for V~i and V^ 1 . Recall that 5, is known because it has been computed 
in the case of order equations. 
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Similarly, equation (j92[) do not involve small divisors nor obstructions. 
However, equation (|91j) does have obstructions and small divisors. In order 
to overcome this problem, we denote by F and G the solutions of 

-F(0)- t iF(0 + oj) = H^(0), 

-G{0) - iiG(6 + u) = D\{0) 
A 1 



D\{9) = -E\{9) + 8^C\{0). 



where 
This gives 

vl(e) = 8F(e) + G{9). 

Taking the average of equation (|9ip . we have that 



D\ + SH$ - B°F8 - B°G - Btyg = 0, (93) 

so we can solve for 8 provided that H§ - B°F ^ 0. 

The other orders do not have any problem. For s n , with n > 2, we have 

n 

V?{0) - fi n V?(6 + uj) + Y,^ n ~ k (0)V 2 k (0) = -E%{9) + 6H?-\0) + 8„C?(0), 

k=0 

V?{6) - n n V?{6 + u) = -E^{6) + 8H%~\0) + 5^(0), 

n 

»V?(9) - fV?{e + u:) + Y, B n - k (O)V 4 k (0) = -E%{9) + 8H^\9) + 8^(0), 

k=0 

-V£(0) - fi n V£(0 + u>) = -E2(0) + 8Hl-\0) + 5^C2(6), 
A 1 

(94) 

and equations (fMj) can be solved uniquely for V™, F 3 " and V^ n , for 
n = 2, . . . , L, where L is the degree for the Taylor expansion. Hence, we 
have obtained 8,8^ E R and 

L 

V{9,s) = Y J V n {Q)s n , 

n=0 

so that the improved solution is 

W <- W + MV, 
A <r- A + 8, 

(JL4- (JL + 6^. 

Remark 18. For L = 1, the algorithm allows us to compute simultaneously 
the invariant torus and the associated linear subspaces. 

The algorithm for the simultaneous computation of the whiskers and the 
invariant torus is 
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Algorithm 6 (Computation of tori and whiskers simultaneously). Con- 
sider given F , ui, Kq and a fixed order L. Given an approximate solution 
(W,[i,fj,), perform the following calculations 

1. Compute E(9,s) = FoW{8, s) -W(9+uj, fis) -((JoK )~ 1 DK )(9 + 
u)n 

2. Compute 

(2.1) a(6,s) = D W(8,s) 

(2.2) 0(9, s) = d s W(9,s) 

(2.3) N(9,s) = [a(9, s) L a(9, s)}- 1 

(2.4) N(9,s) = [(3(9,s) ± p(9,s)]- 1 

(2.5) P(9,s) = a(9,s)N(9,s) 

(2.6) Q(9,s) = (3(9,s)N(9,s) 

(2.7) 1 (9,s) = (JoW(9,s))- 1 P(9,s) 

(2.8) n(9,s) = (Jo W(9,s))~ 1 Q(9,s) 

(2.9) M(9,s) = [a(9,s) \ 7(6,3) \ (3(8, s) \ n(9,s)} 
(2.10) [M(9,s)]- 1 

3. Compute 

(3.1) C(9,s) = M~ 1 (8 + uj,lis)(J(Kq(9 + uj)))- 1 DK (9 + u) 

(3.2) E(9,s) = M-\9 + uj,hs)E(9,s) 

(3.3) H(9, s) = M- 1 (9 + u, fxs)a(9 + oj, (is) 

4. Compute 

(4.1) A(0, s) = P(9, s^KDF o ^)((9, s )j(9, s) - 7 (9 + w, ^5)] 

(4.2) 5(0, s) = Q(0, s)- L [(Z)F o W)(0, ^(0, s) - ^(0 + w, /is)] 

5. (5.1) Solve for 5, satisfying 



<V = 



(5.2) iSofoe /or V 2 satisfying 

V 2 ° - V 2 ° o T w = -£ 2 ° + C 2 % 

Sfei y 2 ° suc ^ ^ e average is 0. 

6. (6.1) Compute A°(9)V 2 °(8) 

(6.2) So/ve /or V^ satisfying 

[ E\- [ 01(8)5,+ [ A°V 2 °+\[ A°]v 2 ° = 

(6.3) Set V 2 ° = V 2 ° + V 2 ° 

(6.4) Solve for V® satisfying 

V?-V?oT u = -E%-A°V? + C$8 lt 

(6.5) Normalize so that J^ e V^ = 

7. SoZue for V® satisfying 

1 y 4 ° - y 4 ° o r w = -El + els, 



/i 
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8. Solve for V 3 satisfying 

- V 3 ° o T w = -E° 3 + C 3 % - B°V 4 ° 

9. (9.1) Solve for F satisfying 

—F - fiF oT w = Hi 

(9.2) Solve for G satisfying 

—G - fiG oT u = -E\ + S^Cl 
A* 

(9.3) Solve for 5 satisfying 

(-7% + s^dj - Wg - B l vf} + s(h$ - WW) = o 

(9.4) Set Vl = 5F + G 
10.(10.1) Solve for V 3 satisfying 

fiVj - fiVi o T u = -El + SH° 3 + 5^1 - B°Vi - B 1 ^ 

(10.2) Normalize so that Jj, e V 3 = 

(10.3) Solve for V 2 satisfying 

Vj - fiVi o T w = -E\ + 5H% + 

(10.4) Solve for V± satisfying 

Vl - fiVi o T u = -E\ + 5Hl + 8pC{ - A°Vi - A*V? 

11. For n = 2 . . . L do 

(11.1) Solve for V 2 n satisfying 

V 2 n - v n V 2 n o T u = -E^{9) + SH?- 1 + <5 M C 2 " 

(11.2) Compute 

n 

A n = s £ j A n - k V 2 k 

k=0 

(11.3) Solve for V™ satisfying 

Vi - n n V? ° T w = -E'{ + 6H?- 1 + SpC? - A n 

(11.4) Solve for V 4 satisfying 

A* 

(11.5) Compute 

n 

B n = B n - k Vf 

k=0 

(11.6) Solve for V 3 satisfying 

»V 3 n - fi n V 3 n o T w = -E$ + 5H%- 1 + SpCS - B n 
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12. Compute 

L 

v{e,s) = Y J v n {e)s n 

n=0 

13. Set W + MV 

(j, <— fj, + 5 
H i- fi + 5^ 

5.3. A Newton method to compute the whiskers. Assuming that we 
have computed an invariant torus K{9) with the associated stable direction 
V s {9) (resp. unstable direction V u (9)) and the rate of contraction (resp. 
expansion) n, we can use a Newton method to compute the whiskers. We 
concentrate on the case of maps, referring to Appendix [E] for the flows. 

We consider the invariance equation (|17p . and we assume that we have 
an initial approximation W for the whiskers, expressed as a power series 

oo 

W(0,s) = Y J W n {9)s n 

n=0 

and such that 

W°(6) = K{9) and W l {9) = V s (9) 
(the unstable case is analogous). 

Remark 19. Again, as discussed in Remark [T5l we do not need to consider 
the modified invariance equation (|18h with the counterterm. The fact that 
K {&) is a solution of equation (jSJ) for exact symplectic maps, implies that 
A = in (I18p . This is guaranteed by the vanishing Lemma in [FdlLS09a]. 

Then, it is clear that the error E for the initial approximation W is such 
that _ 

E{e,s) = Y u E n (0)s n , 

n>2 

since this is exact for the terms of order and 1. 
For a given function G : T £ x R -> M, we denote 

G(9, s) = G [<L] (9, s) + G [ ^ L] (9, s) (95) 

where 

L— 1 oo 

G [<L] (9,s) = Y,G n (9)s n , G^ L] {9,s) = ^ G n {9)s n . 

n=0 n=L 

Using this notation, the linearized equation for the Newton method is 

[DFoW{9,s)]A^(9,s) - ^ 2 \9 + oo^s) = -E^ 2 \9,s). 

Similarly as we did in the previous section, we can perform the change of 
coordinates given by the matrix M(6, s) in ()82[) and reduce the problem to 
solving for V(9, s) the following equation, which is diagonal in Fourier- Taylor 
series, 

R(9, s)V^ 2 \9, s) - V [ ^ 2] (9 + w, (is) = -E [ - 2] (9, s), 
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with R(6,s) given in JBSJ), E(9,s) = M(9 + u, fj,s)- 1 E(9,s) and A = MV. 

Notice that in this case, we do not have to solve the system of equations 
for order and 1 and we can go straight to order n > 2. We use sub indices 
to denote coordinates and superindices to denote the order. Hence, for order 
n > 2 we need to solve the system of equations 

n 

V?{0) - fi n vp(e + u) + Y, An ~ k (0)V 2 k (9) = -E?(9), 

k=2 



(96) 



V 2 n (9)-n n V 2 n (9 + uj) = -E2(9), 

n 

fj,V 3 n (9) - (x n V 3 n (9 + cj) + Y,B n - k (0)V 4 k (0) = -El, 

k=2 

-V£(9)-ix n V?{9 + oj) = -E2. 
A* 

The solution of (|96p for n = 2, 3 provides an exact solution of the invari- 
ance equation up to order 4. That is, if we set 

V l<i] (9,s) = V 2 (9,s) + V 3 (9,s) 

where V 2 and V s are obtained by solving equations (I96p , then the improved 
solution W is given by 

W(9,s) = W(9,s) + M(9,s)V [<4] {9,s), 

where M(9,s) was introduced in (|82[) . The function W approximates the 
solution of the invariance equations with an error E such that 

E(9,s) = E^(9,s). 

This process can be iterated and at each step we solve the invariance 
equation exactly up to an order which is the double of the one we had for 
the initial approximation. Thus, if we assume that the initial guess W is 
such that the error in (1781) satisfies 



E = E^ L \ 

then the modified linearized equation for the Newton method is such that 

R(9,s)V^(9,s) -V^ L \9 + u^s) = -E^ L \9,s), 

with R(9,s) given in (|83j) . If we solve the system of equations (f96|) for 
n = L . . . (2L — 1), then the improved W is 

W{9,s) = W(9,s) + M{9,s)V [<2L] {9,s), 

with M(9, s) as in (|82p . and the new error E satisfies E(9, s) = E^ 2L \9, s). 
The algorithm in this case is: 
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Algorithm 7 (Computation of whiskers). Given F, uj as well as K,V S , 
and an approximate solution W such that 

F{W{9,s)) - W(6 + oj,iis) = E [ ^ L] (9,s) 

with L > 2 and W(9,0) = K(9) and d s W(9,0) = V s (9). Perform 
following calculations: 

1. Compute E^ (9, s) = F o W{9, s) - W{9 + u, fis) 

2. Compute 

(2.1) a(9,s) = D e W(9,s) 

(2.2) (3(9, s) = d s W(9,s) 

(2.3) N(9,s) = [a(9, s) x a(#, s)]- 1 

(2.4) N(9,s) = [(3(9, 3^(3(9, s)}' 1 

(2.5) P(9,s) = a(9,s)N(9,s) 

(2.6) Q(9,s) = p(9,s)N(9,s) 

(2.7) 1 (9,s) = (JoW(9,s))- 1 P(9,s) 

(2.8) V (9,s) = (JoW(9,s))- 1 Q(9,s) 

(2.9) M(9,s) = [a(9,s) \ j(9,s) \ (3(9, s) \ V (9,s)} 
(2.10) [M(9,s)]- 1 

3. Compute 

E^ L \9, s) = M~ 1 (9 + u, [is)E^ L \9, s) 

4. Compute 

(4.1) A(9, s) = P(9, s)^[(DF o W)(9, a) 7 (0, s) - 7 (0 + u, »s)\ 

(4.2) s) = Q(9, s^KDF o ^(0, s) - I^((9 + w, HI 

5. For n = L . . .2L — 1 do 

(5.1) .Sofoe /or V 2 n satisfying 

V 2 n -^ n V 2 n oT UJ = -E^(9) 

(5.2) Compute 



A n = Y^ A n - k V 2 k 

k=L 



(5.3) So/we for V™ satisfying 

V? - [i n V? oT u = —Ei - A n 

(5.4) Solve for V™ satisfying 

-V£ - fi n Vl l o T w = -El 
A* 

(5.5) Compute 

n 

B n = Y^ B n - k Vf 

k=L 

(5.6) Solve for V% satisfying 

»V£ - o T w = -El - & 
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6. Compute 

2L-1 

V(9,s)= V n (9)s n 

n=L 

7. SetW^W + MV 

Appendix A. Fast algorithms to solve difference equations 

WITH NON CONSTANT COEFFICIENTS 

In this section we present fast algorithms to solve for A(9) the cohomology 
equation with non constant coefficients 

A{9)A{9)-A{9 + uj)B{9) = r](9) (97) 

for given A(9), B{9) and rj(9) satisfying either \\A\\ < 1, < 1 or 

HA- 1 ! < 1, < 1. 

Equations of this form appear in the Netwon step for whiskered tori (See 
the informal description in Section [T]). Equations of this form also appear 
in the calcultion of the invariant splitting (see (|67|) - (|68p ). 

We will present two algorithms. The first one is an iterative method with 
an accelerated convergence and the second one very fast (see Section IA.1|) . 
The second one is only for the case of one-dimensional bundles and it is 
faster (computations are 0(N))(see Section [A, 2p . 

A.l. Fast iterative algorithms for the cohomology equation. In this 
section we will present a fast algorithm to solve equation (J97]) using iterative 
methods. We refer the reader to [HdlLS09] where a similar idea is used to 
compute iteration of cocycles. 

We consider first the case < 1 and ||5|| < 1 or, more generally, 

HA -1 ^)!! • \\B(9)\\ < 1 Then, multiplying (JgZj) by A~ 1 {9) on the LHS, we 
obtain 

A(0) = A- 1 (9)A(9 + u)B(9) + A- 1 (9) v (9). (98) 

This is a contraction mapping and it is straightforward to iterate it an obtain 
an algorithm that converges faster than exponentially. 

Next, we compute A(9 + to ) by shifting (|98p and substituting again in 
([98]) . so that we get 

A(0) = A -1 (%(0) 

+ A-\e)A- x {9 + uj) v (9 + uj)B{9) 

+ A~ X {9)A~ 1 {9 + u)A(9 + 2uj)B{9 + lo)B(9). 

Notice that if we define 



fj(9) = A~Ho)v(e) 
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and 

A^{9) = A- 1 {e)A- 1 (9 + 
B 1 {6) = B(9 + u)B(9), 

m (e) = fj(p) + A^(e)f]{e + u)B(o), 

we have that 

A(0) = m(0) + A^(9)A(e + 2oj)B 1 (9) 

which has the same structure as (|98p and we can repeat the same scheme. 
This leads to an iterative procedure to compute A{9), converging superex- 
ponentially in the number of iterations. Thus, define 

A- l +1 {e) = A-\0)A-\6 + 2 n u), 

B n+1 (0) = B n {9 + 2 n u)B n {9), 

Vn+i{9) = Vn(0) + A" 1 (9) Vn (9 + 2 n Lo)B n (9), 

for n > 0, with 

A^(6) = A-\9), 
B (6) = B(8), 
m (9) = fj(9)- 

Then, 

A(0) = Vn+1 (9) + A~l 1 {9)A{9 + 2 n+1 uj)B n+1 (9), V n > 

and 

A(9)= lim r] n+1 (e). 

n— >+oo 

The convergence of the algorithm is guaranteed by the contraction of A' 1 
and B. The cost of computing 2™ terms in the sum is proportional to n 
since it involves only n steps of the algorithm. 
The iterative algorithm is the following: 

Algorithm A. 8 (Solution of difference equations with non constant coef- 
ficient). Given A{6), B(9) such that \\A- X (9)\\ ■ \\B{9)\\ < k < 1, and r/(0) 
perform the following operations: 

1. Compute A(9) = A- 1 {9)r ] {9) 

2. Compute 

(2.1) A(0) = A' 1 {9)A{9 + uj)B{9) + A(0) 

(2.2) A- 1 (6) = A- 1 (9)A' 1 (9 + uj) 

(2.3) B{9) = B(9 + uj)B(9) 

3. Set A -> A 

A^A 
B -> B 
2ui — > oj 

4. Iterate points 2 — 3 
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The case when \\A\\ < 1 and \\B 1 || < 1 can be done similarly. In this 
case, we multiply fl2U) by B~ l (6) on the LHS so that we obtain 

A(0 + u) = A(9)A(9)B- 1 (9) - •q{6)B~ 1 (9). 

Before applying the iterative scheme we shift by —uj. In this way, we have 

A(0) = A{e')h{9')B- 1 {e') - r,{9')B- l (9') 

where 9' = T_ w 0. 

The iterative algorithm in this case is 

Algorithm A.9. Given A(9), B{9) ||^4_(6>) || H^B - 1 (6>) || < k < 1 and r](9), 
perform the following operations: 

1. Compute A(0) = -rj(9)B- 1 (9) 

2. Compute 

(2.1) A(0) = A(0)A(9 - oj)B- 1 (9) + A(0) 

(2.2) A{9) = A(9)A(9 - uj) 

(2.3) B-\9) = B~ 1 {9-uj)B- 1 {9) 

3. Set 

A ->■ A 

S B 

2w — >• a; 

4. Iterate parts 2-3 

This algorithm gives A(0 + uj). Shifting it by — uj we get A(0). 

A. 2. Fast algorithm for the 1-D cohomology equation with non- 
constant coefficients. In this section we present an efficient algorithm for 
the one-dimensional version of equation (|97p . It is an adaptation of methods 
used in [Her83]. 

Consider the following equation, 

which is obtained from (|97p multiplying by B _1 {9) (recall that in this case 
B(0) is just a number). 

We will solve ([99]) in two steps: 

1. Find C(0) and i/£R such that 

AOS) C(9) . . 

im = v WTZT) (100) 

2. Applying (llOOp in (]99p and multiplying by C{9 + uj) we obtain 

z^C(0)A(0) - C(9 + uj)A{9 + uj) = fj(9) (101) 

where fj(6) = C{9 + uj)B- 1 (9)rj(9). 

If we change the unknowns in (llOip by W(^) = C(0)A(0), we are left 
with the equation 

uW(9)-W(9 + uj) = f/(9). (102) 
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Of course, if \v\ 7^ 1, equation (jl02H can be solved in Fourier space. That is, 
we can obtain the Fourier coefficients of W as: 



and the solution is unique. Notice that whenever \u\ = 1, equation (1102p 
involves small divisors, which is not the case for the iterative methods that 
will be discussed in the following section. 
Finally, once we have W{9) we get 

A(0) = C~ l {9)W{9). 

Step 1 can be achieved by taking logarithms of (|1U0|) . Assume that 
A(9)/B(9) > 0, otherwise we change the sign. Then, we get 

log A(9) - log B{9) = logz/ + logC(0) - logC(0 + 

Taking log v to be the average of log A{9) — log -6(0), the problem reduces 
to solve for logC(0) an equation of the form fj23[) . Then C{9) and v can 
be obtained by exponentiation. Notice that logC(0) is determined up to a 
constant. We will fix it by assuming that it has average 0. 

Hence, we have the following algorithm: 

Algorithm A. 10 (Solution of difference equations with non constant coeffi- 
cient (ID)). Given A{9), B{9) andn(9). Perform the following instructions: 

1. (1.1) Compute L{9) = \og{A{9)) - log(£(0)) 
(1.2) Compute L = J* T< L 

2. Solve for Lq satisfying 

L c (9)-L c oT UJ (9) = L(9)-L 

as well as having zero average. 

3. (3.1) Compute C{9) = exp(L c (0)) 
(3.2) Compute v = exp(L) 

4. Compute fj(6) = C(9 + u)B- 1 {9)n{9) 

5. Solve for W satisfying 

vW(9) - W(9 + u) = fj(9) 

6. The solution of (J97D is A(0) = C-\9)W{9) 

Appendix B. Fast Newton method for whiskered tori in 

HAMILTONIAN FLOWS: THE CENTER DIRECTIONS 

In this section, we provide the numerical algorithm to solve the invariance 
equation (fTTI) and the modified one (fT3l) using a Newton method analogous 
to the one described in Section 14.21 

The automatic reducibility can also be proved in this context (see [FdlLS09a]) 
and we provide here the algorithm only. 
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Algorithm B.ll (Newton step for flows in the center direction). Consider 
given X = J(K)VH, oo, Kq and an approximate solution K (resp. K,\). 
Perform the following calculations 

1. (1.1) Compute d w K. 

(1.2) Compute XoK (1.3) Compute the invariant projections n c ,Il M ,H 

2. Set E c = W(d u K -XoK) (resp. set E c = U c {d u K - X o K - ( J o 
K Q )~\DX o K Q )\)) 

3. Following gQJ) 

(3.1) Compute a{8) = DK{6) 

(3.2) Compete /3(0) = J(K (9))- 1 a(9) 

(3.3) Compute 0(9) = a(9)N(8) 

(3.4) Compute 7 (0) = ( J(K(9)))- 1 f3(9) 

(3.5) Compute M{9) = [a(6) \ 7 (0)j 

(3.6) Compute M(9 + w) 

(3.7) Compute (M{9 + u) ± J( y K(9 + uj))M{9 + w))" 1 

(3.8) Compute E{9) = {M{9 + u) ± J{K{8 + w))M(0 + w)) _1 .E(0) c 

(3.9) Compute 



S(9) = P x (6)(ld 2d -m<x(9r){DX{K{e)) + DX(K(9)) )/3(9) 



4. (4.1) Solve for W2 satisfying 




(resp. 



(4.1') SWue /or 0" satisfying 




(4.2') SWue /or W2 satisfying 



E 2 + B 2 5 



Set W 2 such that its average is 0.) 
5. (5.1) Compute S(9)W 2 (9) 



(5.2) .Sofce /or satisfying 



(5.3) Find Wi solving 



[ E l {6)+ f S(9)W 2 (9) + f 



S(0) W 2 = 



d^W x = —Ei - S{W 2 + W 2 ) 



Normalize it so that fj e W\ = 
(resp. 

(5.1') Compute S{9)W 2 {9) 
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(5.2') Solve for W 2 satisfying 



E x {6) + I B x {9)8- I S(6)W 2 {0) 



S(0) 

Jt 



W 2 = 



(5.3') Find W\ solving 

d u Wi = -Ei - S(W 2 + W 2 ) + B16 
Normalize it so that Jj t W\ = 0). 

6. The improved K is K(0) + M(9)W(6) 
(resp. the improved X is A + 5). 

Notice t hat st eps (1.1), (3.1), (4.1) (resp. (4.2')), (5.3) (resp. (5.3')) in 
Algorithm IB. Ill are diagonal in Fourier series, whereas the other steps are 
diagonal in the real space representation. The algorithm only stores vectors 
which are of length of order N. 

Appendix C. The order by order method for whiskers in 

Hamiltonian flows 

In this section we present the result analogous to the one described in 
Section 15. II to solve the invariance equation (|17p for the whiskers in the case 
of Hamiltonian flows. 

As in Section 15.11 we look for W as a power series 



W(9,s) = J2w n (9)s n , 



n=0 

and match similar coefficients in s n on both sides of equation (|17p . 
For n = 0, one obtains 

d u W (9) = {XoWo)(9) (103) 

which admits the solution Wq(8) = K(9), where K is a parametrization of 
the invariant torus. 
For n = 1, we obtain 

SL,Wi(0) + Wi(0)/z = (DX o K(fi))Wi(e), (104) 

from where we read that W\{6) is an eigenfunction with eigenvalue — /i of 
the operator C u 

Cw :=d w -DXoK{e). 
Again, we note that, multiplying a solution of (|104p by a scalar b E R, we 
also obtain a solution. See Remark 1121 
For n > 2, we obtain 

d u W n (0) + W n (6)nfi = (DX o K(8))W n (9) + R n (W , . . . , W^), (105) 

where R n is an explicit polynomial in Wq, . . . , W n -\ whose coefficients are 
derivatives of X evaluated at Wq = K. 



(106) 
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Notice that, in this case, equation (|105|) can be solved provided that n/i 
is not in the spectrum of the operator (this is a non-resonance condition 
which is clearly satisfied since the stable spaces are 1-dimensional) . As in 
the case of maps, the previous equation can be solved using the large matrix 
method. 

Appendix D. A Newton method to compute simultaneously the 

INVARIANT TORUS AND THE WHISKERS FOR FLOWS 

As in Section 15.21 we start with an initial approximation (W, fi) (resp. 
(W, //, //) for the invariance equation (fT7|l (resp. (fTTl) ). that is 

X(W(6, s)) - (d u + pus^A W(9, s) = E(9, s), 

X(W(6, s)) - (d u + ns-^j W(9, s) - (J~ X DX) o KqX = E{6, s) 

and we will look for an improved solution 

W -> W + A 
A -> A + <5 

fJL -t /i + <5 M 
by solving the linearized equation 

{DX o (W(9, s))A(6, s) - (d w + fMs^j A(9, s) - s-^-W(9, s)6 

-(J- 1 DX)oK 5^ = -E(6,s). 

Once again, we will use a reducibility argument similar to the automatic 
reducibility of Lagrangian tori. This will lead to a diagonal equation. Ap- 
plying the operators Dq and d s to equations (|lU6p . we obtain 

DX(W(9, s))D e W(9, s) - (d w + fis^j D e W(9, s) = 0(E), 

DX(W(9, s))d s W(9, s) - (d w + fis-^-j d s W{9, s) = fid s W(9, s) + O(E) 

The vectors ( J o W^DgWN and (J oW)~ l d s W N , where N and N are 
normalization matrices (see (I8ip ). are the symplectic conjugate vectors of 
DgW and d s W, respectively. 

By the Hamiltonian character of the vector field, we have that 

(DX o W(9, s))((J o W)- 1 D e WN)(9, s) - (d u + fis^j ((J o W^DgWN)^, 

S(9,s)D e W(9,s) + 0(E) 
(DX o W(9, s))((J o W)- 1 d s WN)(9, s) - (d u + ns-^\ ((J o W)- 1 d s WN)(9, s 

- /x((J o W)- 1 d s WN)(9, s) + B(9, s)d s W(9, s) + 0(E) 
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where S(9, s) and B(9, s) are matrices which can be computed. 
We summarize these properties in the following proposition: 

Proposition D.12. Using the same notations (18ip as in Proposition^ and 
considering the matrix M(0,s) introduced in (|82p . we have 



DX o W(9, s)M(e, s) 
where 

11(6, s) 

with 



duj + fJ-s^ ) :U(0..s) - M (f).. s)K(f)..s) + 0(E) 



S{9,s) 



o 


O 


fi B(6,s) 
-/x 



(108) 



S(0, s) = P L (6, s)[d u >y(e, s) - DX(W(6, s)) 7 (9, s)] 
B(0, s) = Q ± (0, s)[d u ri(6, s) - DX(W(6, s))n(6, s)] 

or, equivalently, 

S(8) = P ± (9, s)(Id -P(6, s)a{9, s) L )(DX(W{0, a)) + DX{W{9, s)) ± )P(0, s) 
B(9) = Q x (9, s)(Id -Q(9, B )P(0, s) ± )(DX(W(9, s)) + DX(W(9, s)) ± )Q(0, s) 

and E is the error in (|78p . 

Now if we change the unknowns A = MV in (|107j) and multiply by 
M~ 1 (e,s) the LHS, by Proposition EEJ we are left with the following 
system of equations 



K{e,s)v{e,s) 



+ Ms^J V(0, s) - C(9, s)6 u = -E(9, s) + s5H(8, s) 

(109) 



where lZ(9,s) is given in (|108p and 

C(0,s) = M- 1 (9,s){J- 1 DX)oK {9), 
E(9,s) = M- 1 (8,s)E(6,s), 
H{9, s) = M~ l (9, s)d s W(8, s). 



We expand the terms in (|109p as a power series up to some order L and 
match coefficients of the same order on both sides of the equation. We 
use subindices to denote coordinates and superindices to denote the order. 
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Hence, for order s° we have 

- d u V?(0) + s°(e)v 2 °(e) - Cl{9)5^ = -E%(6), (110) 

- d u V§{9) - Cl(9)5^ = -E° 2 (8), (111) 

/^ 3 » - d„v*(e) + B\e)v2(e) - c!(9)5, = -e° 3 (9), (112) 

- vVgip) - 8^(8) - C°(0)<^ = -E° A (8). (113) 

Notice that (jllOp and (jllip can be solved using the Algorithm IB. 1 1 1 
Hence, we determine V®, V 2 and 5^. Once we know 6^, we can solve uniquely 
for V 3 ° and V® in equations (|112p and (| 1 1 3 [) . These equations do not have 
any small divisors nor obstructions. 

For order s 1 we have 

- d u V?(0) - uV{(9) + S\9)V 2 \9) + S\9)V 2 °(9) (114) 

= -El(8) + 5H°(9) + 5^1(8), 

- d u Vi(0) - nV 2 \9) = -E l 2 {9) + 5H$(6) + 5^(9), (115) 

- d u V£(0) + B°(9)Vi(9) + B\9)Vi(9) (116) 

= -El(8) + 5Hl{9) + 5,Cl(8), 

- d u Vl(6) - 2fiVl(6) = -E\{8) + 8Hl(9) + 8^C\{9). (117) 

Notice that once we choose 5, equations (j!14j) and (|115p are uniquely 
solvable for V/ and V 2 . Recall that 5^ is known, since it has been computed 
in the case of order equations. 

Similarly, equation (|1 17[) can be solved without small divisors nor obstruc- 
tions. However, equation (|1 16[) does have obstructions and small divisors. 
In order to overcome this problem, we denote by F and G the solutions of 

-d u F(0) - 2nF{6) = Hl{9), 
-d U3 G(9)-2fiG(9) = Dl(9) 

where 

Dl(9) = -E\{8) + 5^Cl(9), 

then 

V A 1 {9) = 5F{9) + G{8). 
Taking averages of the equation for V^ 1 we get 

Dj + SHji - WFS - WG - B l Vg = 0. 

So we can solve for 5 provided that — B°F ^ 0. 
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Now the other orders do not have any obstructions, 

n 

- d u V?(0) - nnV?{6) + S n ~ k {9)V 2 k {9) = -E?(6) + SH^-\e) + 6^(9), 

- d u V?(0) - n^V 2 n (9) = -E%{9) + 8H^ 1 (9) + 5^(9), 

n 

- d„V£{9) - (n - l)pV?{0) + £ B n -\9)V k {9) = -E%{9) + 8H^- 1 (9) + 6^(9), 



k=0 



- d u V?(9) - (n + 1)^(9) = -El(9) + 5H2~ L (9) + 5^2(9). 

(118) 

for n > 2 and they can be solved uniquely for V™, V 2 n , V£ and V£", for 
n = 2, . . . , L, where L is the degree for the Taylor expansion. Hence, we 
have obtained 5, 5^ € R and 

L 

V(9, S ) = Y,V n (9)s n 



n=0 



and the improved solution is 



W <- W + MV 
fi <- n + 5^ 

The algorithm for the whiskers and the invariant torus, analoguous to 
Algorithm [6l is 

Algorithm D.13 (Computation of whiskers and tori for flows). Consider 
given X , u, Kq and a fixed order L. Given an approximate solution (W, /x, fi), 
perform the following calculations 

1. Compute E{9,s) = X(W(9,s)) - (d u + fisd s )W(9,s) - (J^DX) o 

2. Compute 

(2.1) a(9,s) = D e W(9,s) 

(2.2) P(9,s) = d s W(9,s) 

(2.3) N(9,s) = [a(9, s) x a(#, s)}- 1 

(2.4) N(9, S ) = if3(9,s)^P(9,sT 1 

(2.5) P(9,s) =a(9,s)N(9,s) 

(2.6) Q{9,s) = (3{9,s)N{9,s) 

(2.7) - / (9,s) = (JoW{9,s))- 1 P(9,s) 

(2.8) rj(9, s) = (J o W(9, s))~ 1 Q(9, s) 

(2.9) M(9,s) = [a(9,s) \ j(9,s) \ 0(6, s) \ r](9,s)\ 
(2.10) [M(9,s)]- 1 

3. Compute 

(3.1) C(9, s) = M~ 1 (9, s)(J- x DX) o K (9) 



52 G. HUGUET, R. DE LA LLAVE, AND Y. SIRE 

(3.2) E{9,s) = M- 1 (8,s)E(6,s) 

(3.3) H(6,s) = M- 1 (e,s)(3(9,s) 

4. Compute 

(4.1) S{9, s) = P ± {9, s)(Id -P{9, s)a(9, s) ± )(DX{W{9, s))+DX(W{9, s)) ± )P(9, s) 

(4.2) B(9, s) = Q ± (9, s)(Id -Q(9, s)p(0, s) ± )(DX(W(9, s))+DX(W{9, s)) ± )Q(0, s) 

5. (5.1) Solve for 5^ satisfying 

! E° 2 -\ [ C 2 °l S„ = 
(5.2) So/ue /or V 2 satisfying 

-d u v 2 ° = -El + c 2 % 

5ei suc/i that the average is 0. 

6. (6.1) Compute S°(9)V 2 °{9) 

(6.2) Sofee for V 2 satisfying 

[ E°- [ C»{9)5,+ [ S°V 2 °+\[ S°\v 2 ° = 

(6.3) Set V 2 ° = V 2 ° + V 2 ° 

(6.4) Solve for V® satisfying 

-d„V? = —Ei - S°V 2 ° + C% 

(6.5) Normalize so that fj, e V^ = 

7. Solve for V® satisfying 

-A*V? - 9 w Vf = + C26„ 

8. SoZue /or Kj° satisfying 

VV 3 ° - d u V$ = -El + C 3 % - B°V2 

9. (9.1) Solve for F satisfying 

-d u F - = Hi 

(9.2) Solve for G satisfying 

-8 U G - 2^G = -El + S^Cl 

(9.3) Solve for 5 satisfying 

[~E\ + 8^dJ - WG - B l vf} + 5(H$ - WF) = 

(9.4) Set V} = 5F + G 
10.(10.1) Solve for V% satisfying 

-d^Vi = -E\ + <5#3° + 5^Cl - B°Vi - B x v2 

(10.2) Normalize so that f JeU V% = 

(10.3) Solve for V 2 satisfying 

-d u Vi - fiVi = -E\ + 5H Q 2 + b^Cl 
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(10.4) Solve for V± satisfying 

—d u V-l - iiVl = -E\ + 5H° X + 5^C\ - S°Vi - S 1 ^ 

11. For n = 2...L do 

1) Solve for V£ satisfying 

-d u V? - niiV? = -E%{9) + SH?- 1 + 
1) Compute 

n 

s n = Y^ S n ' k V 2 k 

3) Solve for V™ satisfying 

-d^ - nnV? = -E% + 5H?- 1 + S^Ct - S n 
1) Solve for V£ satisfying 

-d u V? - (n + l)/xy 4 n = -El + SH2' 1 + b^Cl 
5) Compute 

n 

3) Solve for V£ satisfying 
-d„V£ - (n - l)fiV 3 " = -El + 5H™- 1 + 5^C% - B n 

12. Compute 

L 

V{9) = Y^V n (9)s n 

n=0 

13. Set W^W + MV 

A <- X + 5 
fi <- fi + 5^ 

Appendix E. A Newton method to compute the whiskers for 

flows 



(11.1) 
(11.2) 

(11.3) 
(11.4) 
(11.5) 

(11.6) 



We consider the invariance equation (|T7|) for flows and we assume that 
we have an initial approximation W for the whiskers, expressed as a power 
series 



W(9,s) = J2w n (9)s r - 



n=0 

and such that 

W°(9) = K(9) and W 1 (9) = V s {9) 

(the case unstable is analogous). 

Then, it is clear that the error E for the initial approximation W is such 
that _ 

E{9,s) = Y J E n {9)s n 

n>2 
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because the approximation is exact for the terms of order and 1. 

Using the notation introduced in (195p . the linearized equation for the 
Newton method is 

[DX o W(0, s)]A^ 2 \0, s) - (d u + »sds)A^(0, s) = -E^ 2 \d, s). 

Proceeding as in the previous section we can perform the change of coordi- 
nates given by the matrix M(9, s) in (|82p and reduce the problem to solving 
for V(9, s) the following equation, which is diagonal in Fourier- Taylor series, 

11(0, s)V [ - 2] (9, s) - (d u + fisds)V [ - 2] (9, s) = -E [ ^ 2] (9, s), 

with K(9, s) given in (FTPS]) and E(0, s) = M(6, s)- l E(0, s). 

Notice that in this case, we do not have to solve the system of equations 
for order and 1 and we can go straight to order n > 2. We use subindices 
to denote coordinates and superindices to denote the order. Hence, for order 
n > 2, we need to solve the system of equations 

n 

- W W - nfiVr(8) + ^ S n - k (9)V 2 k (6) = -E?(6), 

k=2 

(119) 

n 

- d u v?(o) - (n - i)W(o) + J2B n - k (d)v 4 k (e) = -Eg, 

k=2 

- d u V£(0) - (n + l)fiV?(9) = -El 

Notice that now the solution of (|119p for n = 2, 3 provides an exact 
solution of the invariance equation up to order 4. That is, if we set 

V [<4] (9,s) = V 2 {9,s) + V 3 (9,s) 

where V 2 and V s are obtained by solving equations (|119p , then the improved 
solution W given by 

W(8, s) = W(9, s) + M{9, s)V [<4] (8, s), 

where M(9,s) was introduced in (|82p . satisfies that it approximates the 
solution of the invariance equations with an error E such that 

E(8,s) = E^{9,s). 

This process can be iterated and at each step we solve the invariance 
equation exactly up to an order which is the double of the one we had for 
the initial approximation. Thus, if we assume that the initial guess W is 
such that the error in (I106P satisfies that 

E = E^ L \ 

then the modified linearized equation for the Newton method is such that 
K{0,s)V^ L \9,s) - (d^ + fisds^^iO^) = -E^(9,s), 
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with 1Z(9, s) given in (1108H . If we solve the system of equations (|119p for 
n = L . . . (2L — 1) then the improved W is 

W(9, s) = W(9, s) + M{9, s)V [<2L] (6, s), 

with M{9, s) as in (|82p . and the new error E satisfies E(9, s) = E^- 2L \9, s). 
The algorithm in this case is 

Algorithm E.14 (Computation of whiskers for vector- fields) . Given X, uj 
as well as K,V s ,fJ, and an approximate solution W such that 

X o W(9, s) - (d u + ^sd s ) W(9, s) = E^ (9, s) 

with L>2 and W(9, 0) = K(0) and d s W(9, 0) = V s {6), perform the follow- 
ing calculations: 

1. Compute E^ (9, s) = X o W(9, s) - {d u + fisd s )W(9, s) 

2. Compute 



(2.1) 


a(0, s) 


= D e W(9,s) 


(2.2) 


/3(9, s) 


= d s W(9,s) 


(2.3) 


N(9,s) 


= [a(9, s) ± a(9, s)}' 1 


(2.4) 


N(9,s) 


= [p(9, 3)^/3(9, s)}- 1 


(2.5) 


P(9,s) 


= a(9,s)N(9,s) 


(2.6) 


Q(9,s) 


= /3(9,s)N(9,s) 


(2.7) 


1(0, s) 


= (JoW(9,s))- 1 P(9,s) 


(2.8) 


rj(9,s) ■- 


= (Jo W(9,s))~ 1 Q(9, s) 


(2.9) 


M{9, s) 


= [a(9,s)\ 1 (9,s)\f3(9,s)\ V (9, s 


(2.10) 


[M(9,s 


r 1 



3. Compute 

E^ L \9,s) = M~ 1 (9,s)E^ L \9,s) 

4. Compute 

(4.1) S(9, s) = P x (9, s)(Id -P(9, s)a(9, s) ± )(DX(W(9, s))+DX(W(9, s)) ± )P(9, s) 

(4.2) B(9, s) = Q ± (9, a) (Id -Q(9, s)(3(9, s) ± )(DX(W(9, s))+DX(W(9, s)) ± )Q(9, s) 

5. For n = L . . . 2L - 1 do 

(5.1) Solve for V 2 n satisfying 

-d U] V 2 n (9)-n^V 2 n (9) = -E^{9) 

(5.2) Compute 



s n = Y, s n ~ k v 2 k 

k=L 

(5.3) Solve for V{ 1 satisfying 

- W(#) - n^V 1 n (9) = -El - S n 

(5.4) Solve for VI 1 satisfying 

-d u V?(6) ~(n + 1)^(9) = -El 
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(5.5) Compute 

n 
k=L 

(5.6) Solve for V£ satisfying 

-d u V?(0) - (n - 1)»V?(0) = -El - B r 

6. Compute 

2L-1 

V(9,s) = V n (6)s n 

7. SetW <-W + MV 



n=L 
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